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I. 


Introduct Ion 


SOLVER  is  a  root  solver  that  provides  automatic  compilation, 
loading,  and  execution.  It  uses  a  simple  version  of  Muller's  method 
(see  Ref  .  -44-.  SOLVER  is  intended  for  the  occasional  user  uho  needs 
root';  but  balKs  at  learning  to  use  sophisticated  but  complicated 
solvers,  uhich  usually  also  require  .considerable  additional 
prog-amming .  The  user  need  only  supply  the  function  to  be  set  to 
zero  and  the  number  of  solutions  desired.  SOLVER  is  available  on 
the  wRAY- 1  computer,  and  can  be  obtained  by  typing: 
rfilem  read  1222  vcray  solver (LF)end  /tv 


This  report  corresponds  to  the  SOLVER  version  cf  August  31. 
1979.  Later  versions  of  SOLVER  will  be  stored  in  FILEM  directory 
•Cray  of  user  number  1222.  The  user  should  periodically  checK  the 
date  of  this  file  (filom  hou  1222  .Cray  solver)  to  see  if  the 
program  has  been  updated.  The  file  SOLVER  is  a  LIB  file;  it 
contains  the  latest  sources  as  well  as  the  latest  documentation. 

To  obtain  the  documentation,  type: 

lib  solver(LF)x  sol v/doc(LF)end  /  t  v 
net out  Cuscil  solv/doc  box  nnn  solver  /  tv 


A  more  robust  but  more  complicated  Muller  subroutine  is  MRAF. 
CACM  algorithm  196.  Changes  suggested  by  J.  Traub  and  the 
certification  appeared  in  CACM  January  1968  p.12.  A  FORTRAN 
translation  uas  made  and  some  remaining  glitches  uere  removed  and 
new  leatures  added  by  A.B.  Langdon.  DC  BerKeley  Febrary  1968.  This 
subroutine  is  used  in  the  ROOTS  program  (see  Ref.  2.3).  If  the  demand 
is  sufficient,  this  latter  subroutine  may  be  made  available  as  an 
option  at  some  future  date. 


Other  root  solvers  are  available  on  the  MFECC  CDC-7600:  these 
might  be  transferred  to  the  CRAY  by  a  user.  They  include  ABEROOT.  by 
C.E.  Seyler,  uhich  employs  a  global  Neuton  method,  and  routines 
RPZERO.  CPZERO  for  finding  roots  of  real  and  complex  polynomials, 
developed  at  LASL  (see  Ref.  5). 


II. 


Elementary  SOLVER 


To  run  SOLVER,  the  user  enters  "solver"  'allowed  by  some 
commands  from  the  terminal.  In  this  section.  some  basic  commands 
are  introduced: 

<integer>  - 

Specify  the  number  of  roots  desired  (defaults  to  1  and  cannot 
exceed  1024;  this  is  referred  to  as  n  in  the  rest  of  this 
repor t ) . 

f  [file  name]  - 

Enter  the  function  to  be  solved.  If  the  function  is  to  be 
entered  from  the  terminal,  the  user  should  uait  for  the 
prompt  to  appear.  The  name  of  the  function  is  "f“  and 

the  argument  is  "z".  The  input  format  is  assumed  to  be 
that  of  FORTRAN  statements.  In  the  tty-input  mode,  input 
is  terminated  when  “  f 11  is  the  first  non-blank  character 
followed  by  "■="  or  a  blank,  or  a  slash  is  the  first 
non-blanK  character  of  an  input  line. 

go  Cfile  name  or  11  tty"]  - 

Compile,  load  and  execute.  At  execution  lime,  a  number  of 
parameters  can  be  set  to  override  the  default  values. 

If  tty  is  specified,  the  user  should  uait  for  the  prompt 
to  appear.  Inputs  are  in  the  format  of  a  NAMELIST 
(terminated  by  $) . 

Some  of  these  parameters  are: 

n  -  Number  of  roots  desired  (defaults  to  <integer>). 

max i l  -  Maximum  number  of  iterations  per  root  allowed 
(defaults  to  100). 

Known  -  Number  of  Known  roots  (defaults  to  0). 
h  -  Initial  distance  between  estimated  roots 

(defaul ts  to  0.5) . 

epl  -  Relative  error  tolerance  on  z(i)  (defaults  to 

l.e-12  which  is  also  the  minimum  one  can  specify). 
ep2  -  Error  tolerance  on  f(z(i))  (defaults  to  l.e-20 
which  is  also  the  minimum  one  can  specify), 
of  =  0  =>  Output  results  to  tty  only  (default). 

1  =>  Output  results  to  tty  and  a  disk  file. 

2  »>  Output  results  to  a  disk  file  only. 

3  =>  Do  not  output  results. 

See  Section  III  for  the  name  of  disk  file, 
y  -  Array  for  initial  estimated  values  and  result  of 

roots  (defaults  to  0). 

dbug  =  0  »>  Drop  file  deleted  after  execution  (default). 

1  =>  Drop  file  retained  after  execution. 

xeq  Cfile  name  or  "tty"]  - 

Execute  the  generated  program  again  without  recompiling  it. 

See  the  command  "go"  for  details. 

savef  <file  name)  - 


-2- 


Save  the  function  into  the  file  <file  name)  which  can  be  read 
by  using  the  command  "f"  later. 

Kb  in  - 

Keep  the  binary  1  ibrary  "solv/b/1 "  in  the  active  file  area. 

If  SOLVER  is  being  using  often,  this  step  conserves  computer 
effort  (the  cost  of  extracting  “solv/b/1  11  from  the  1  ibrary 
file  each  time  is  eliminated). 

end  - 

Terminate  SOLVER. 


Example  II. 1: 

Most  common  syntax 
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rout 

rout 

rout 

rout 

rout 

rout 

rout 

rout 

rout 

rout 

rout 

rout 
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ine 

ine 
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ine/user 

ine/user 
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ine 

ine 

ine 

ine 
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ine/user 
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solver  /  2  2 

**-li  root  solver  79.08  *f<* 

extracting  solv/b/l 
c  07/20/79  12:39:32  001222 
>3 
>f 


>go 
comp  i 
FT004 
FT001 

t. 


execu 
root  s 

( 

( 

( 

>end 

al  1  done 


f  =  z**3  -  8. 

1 ing  and  1 oad ing 

-  CFT  VERSION  -  07/14/79  SCHEDULER 

-  COMPILE  TIME  *  0.0065  SECONDS 

ray  loader  version  -  cl21  07/05/79 
l  ing 


2 . 00000000000e+00 
-  1 .000000000000+00 
-  1 . 00000000000e+00 


0 . 00000000000e~0 1  ) 
1 .73205080757e+00  ) 
-1 .73205080757e+00  ) 


Example  11.2: 

All  commands  on  execute  line 


user 

rout ine 
rout  ire 
rout ine 
rout ine/user 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 


solver  20  f  go  end  / 

***  root  solver  79.08 
extracting  solv/b/l 
c  07/20/79  12:39:32  001222 
:  f  *  csm  (  z  ) 

compiling  and  loading 
FT0O4  -  CFT  VERSION  - 
FT0O1  -  COMPILE  TIME  * 

*♦*  crag  loader  version 
exocut ing 
roo‘ s : 

(  0.000000000000-01 
(  -3.14159265359e+00 

(  3.14l59265359e+00 


2  3 
*** 


07/14/79 
0.0061 
-  c  121 


SCHEDULER 

SECONDS 

07/05/79 


0 . 00000030000e-0 1  ) 
9 . 33326735795e-2 1  ) 
-4.28706032686O-36  ) 


-3- 


rout ine 

( 

-6 . 283 1 6530? 1 8e+00 

1 . 557863? 1627e-25 

) 

rout ine 

( 

b.28318530?18e+0O 

-1.815153841 10e- 1 7 

) 

rout ine 

( 

-1.256637061440+01 

6. 9447073 1794e-23 

) 

rout ine 

< 

1 .25663706 144e +01 

-4.06743060262e-20 

) 

rout ine 

( 

-9.4247??96077e+00 

3.36131 166524e-2? 

) 

rout ine 

( 

9 . 4247??960?7e+00 

4. 33 163059989e- 19 

) 

rout ine 

( 

-2 . 5 132?41228?e+0 1 

2. 4470828 1560e-24 

) 

rout ine 

< 

2.5132?41228?e+01 

2 . 468492675?6e-  19 

) 

rout ine 

( 

-2. 1991 148575 le+01 

2. 3 1592585436e-23 

) 

rout ine 

( 

2.  1991 1485?51e+01 

3 . 5416 1507604e-  17 

) 

rout ine 

( 

-3 . 45575 19 1895e+0 1 

-9 . 2 1324326052e-27 

) 

rout ine 

( 

3 . 45575 19 1895e+0 1 

7. 4254485260  le- 19 

) 

rout ine 

( 

-1. 884955592 15e+01 

-4. 6053?034333e-22 

) 

rout ine 

( 

3. 14159265359e+0 1 

-3 . 344892 1670 le-  16 

) 

rout ine 

( 

1.88495559215e+01 

-1. 30075028 143e-26 

) 

rout  ine 

( 

-4. ? 1238898038e+0 1 

1.40250535984e-15 

) 

rout ine 

( 

-4. 398229? 1503e+0 1 

5.5025?927188e-20 

) 

rout ine 

al  1 

done 

If  the  user  wants  roots  in  some  kind  of  order,  and  hopefully 
with  none  skipped,  the  user  must  provide  proper  initial  guesses  for 
the  y  values.  It  is  also  possible  to  ask  the  roots  to  be  reordered 
by  the  user  subroutine  (see  next  section). 


III. 


Intermediate  SOLVER 


It  is  sometimes  desirable  to  find  different  sets  of  roots  of 
a  function  uith  parameters  assuming  different  values.  SOLVER  provides 
this  capability  by  allowing  the  user  to  supply  a  subroutine.  The 
structure  of  the  main  program  consists  of  the  following: 

40  cont inue 

call  <the  subroutine'/ 

if  (  1(30). It. 0  )  call  exit  (  dbug  ) 

call  <mullor  to  find  roots> 

write  <resul ts> 

go  to  40 

There  are  variables  which  are  common  to  the  main  program,  the 
function  and  the  user-supplied  subroutine  for  the  user  to  use) 
tney  are: 

integer  1 (30) 

complex  c(30).  p( 1024. 100) 

These  variables  can  be  read  in  either  from  the  terminal  or  from 
the  input  file  during  execution.  As  the  user  may  have  noticed. 

1(30)  is  the  control  variable  for  looping  to  find  roots.  The 
loop  stops  when  1(30)  <  0. 


It  is  possible  to  "follow"  a  root  as  parameters 
change  by  careful  programming  of  the  option  subroutine,  but  the 
user  must  implement  this  him/herself. 


Several  I/O  units  are  also  available  in  the  subroutine: 

59  -  Terminal  (tty). 

05  -  Input  (exists  only  if  "go  <file  name>"  or  "xeq  <file  name>" 
is  given  -  logical  unit  5  is  connected  to  the  file 
<  f i 1 e  name)) . 

06  -  Output  (file  name  is  “rout"  with  the  current  suffix  and 
a  letter  appended.  Exists  only  if  of=l  or  2). 


Libraries  such  as  DISSPLA.  TV80LIB.  FORTLIB  and  BASELIB  are 
1 inKed  at  load  time  automatically.  However,  the  user  may  call 
subroutines  from  other  libraries;  SOLVER  allows  the  user  to  add 
and/or  delete  libraries  at  load  lime. 


The  following  are  additional  commands  available  in  SOLVER: 
sub  Cfilo  name 1  - 

Enter  the  subroutine.  All  input  formal  and  termination 
rules  are  similar  to  the  command  11  f"  in  Section  II. 
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saves  <file  name>  - 

Save  the  subroutine  into  the  file  <  f  i  1 e  name)  which  can  be 
read  by  using  the  command  "sub"  later. 

al ib  Cl ib  name  . . . ]  - 

Add  libraries  at  load  time.  This  command  is  terminated  by 
the  carriage  return  Key.  the  line  feed  Key.  or  the  escape 
Key . 

dl ib  C 1 ib  name  . . . ]  - 

Delete  libraries.  The  syntax  of  this  command  is  the  3ame  as 
that  of  the  command  “alib". 

<else>  - 

If  the  input  is  the  name  of  an  executable  file.  SOLVER  will 
run  it  as  a  control  lee.  This  is  useful  when  using 
a  text  editor  to  modify  the  function,  the  subroutine  or 
the  input  file,  or  using  DDT  to  debug  the  program. 

In  this  case,  the  rest  of  the  input  1 ine  is  passed  to  the 
controllee.  i.e.  any  command  following  will  be  lost  to 
SOLVER. 


Example  1 1 1 .  1 : 

user:  solver  /  2  2 

routine  :  ***  root  solver  79.08  *** 

routine  :  extracting  solv/b/1 

routine  :  c  07/20/79  12:39:32  001222 

rout ine/user :  >3  sub 

rout ine/user :  :  data  1(30)  /3/.  c(l)  /64./ 

rout ine/user :  :c 

routine/user:  :  1(30)  =  1(30)  -  1 

rout ine/user :  :  if  (  l(30).lt.0  )  return 

routine/user:  :  c(l)=c(l)/2. 

rout ine/user :  :  wr i te(59. 100)  c(l) 

rout ine/user :  :  100  format  (  “c(l)  =  ".  2fl4.6  ) 

rout ine/user :  : / 

rout ine/user :  >f  go  end 

rout ine/user :  :  f  =  z**3  -  c(l) 

routine  :  compiling  and  loading 

routine  :  FT004  -  CFT  VERSION  -  07/14/79  SCHEDULER 

routine  :  FT001  -  COMPILE  TIME  *  0.0136  SECONDS 

routine  :  ***  Cray  loader  version  -  cl21  07/05/79 

routine  :  executing 

routine  :  c(l)  =  32.000000  0.000000 

routine  :  roots: 

routine  :  (  3 . 1 74802 10394e+00  .  0 . 00000000000e-0 1  ) 

routine  :  (  -  1 . 58740 105 197e+00  .  2 . 74945927400e+00  ) 

routine  :  (  -1 .58740 1051 97e+00  .  -2.74945927400e+00  ) 

routine  :  c(l)  ■  16.000000  0.000000 

rout ine  :  roots: 

routine  :  (  2 . 5 1984209979e+00  .  -4. 17757316097e-28  ) 

routine  :  (  - 1 .25992 104989e+00  ,  2. 18224727 194e+00  ) 
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roul ine 
rout ine 
rout  ins 
rout ine 
rout ine 
rout ine 
rout ine 


(  -1 . 25992 104989e+00 

c ( 1 )  =  8.000000 

roots : 

(  2.00000000000e+00 

(  -1.000000000006+00 
(  -1.000000000006+00 
all  done 


-2. 18224727 194e+00  ) 
0.000000 

2 . 44343 158944e-46  ) 
1 .73205080757e+00  ) 
-1 .73205080757e+00  ) 


Example  111.2: 

Deleting  the  default  graphic  libraries  which  are  not  used 
because  no  plots  are  made. 


user 

rout ine 
rout ine 
rout ine 
rout ine/user 
rout ine/user 
roul ine/user 
rout ine 
rout ine 
roul ine 
roul  ine 
rout ine 
rout ine 
rout  ine 
rout ine 
rout ine 
rout ine 


solver  /  2  1 

***  root  solver  79.08  *** 

extracting  solv/b/1 
c  07/20/79  12:39:32  001222 
>dl  ib  disspla  tv801 ib 
>2  f  go  end 
:  f  =  z**2  -  16. 

compiling  and  loading 

FT004  -  CFT  VERSION  -  07/14/79  SCHEDULER 

FT001  -  COMPILE  TIME  =  0.0063  SECONDS 

crag  loadeh  version  -  cl21  07/05/79 
***uarn ing-unsat isf ied  external s*** 
execul  ing 
roots: 

(  -4.00000000000e+00  ,  -2.01948391737e-28  ) 
(  4.00000000000e+00  ,  6. 457183323ble-42  ) 
al 1  done 
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IV. 


Advanced  SOLVER 


This  section  introduces  the  plotting  capability  of  SOLVER. 

A.  Plotting  Roots  vs.  One  Parameter 

SOLVER  can  plot  the  roots  of  a  function  vs.  a  user  defined 
parameter  whose  values  are  to  be  stored  in  x.  a  real  array  of  size 
1024.  To  do  so,  the  user  has  to  set  the  variable  "ploimode"  to  1, 

2.  or  3  depending  upon  the  destination  desired  for  the  plot  (see 
below)  and  must  assign  values  to  x  in  the  subroutine  (see  below). 
Several  parameters  are  included  for  the  user  to  set  at  execution 
time;  they  are: 

plotmode  -  »  0  =>  No  plot  (default). 

1  =>  Fr80  file. 

2  =>  Printer  file  (always  called  "disoul"). 

3  =>  Tektronix. 

If  the  value  read  for  plotmode  is  positive,  the  user 
can  vary  the  value  of  plotmode  between  0  and  a 
positive  intpger  dynamically.  The  user  assigns  a 
positive  number  to  plotmode  only  for  those  roots 
he/she  wants  to  plot.  The  number  of  abscissa  values 
plotted  is  the  number  of  times  the  plotmode  is  found 
to  be  positive  upon  return  from  the  subroutine. 

plot  id  -  Depending  on  the  destination,  plot  id  has  different 
interpretat ions.  If  plotmode  is 

1  (fr30  fiie)  =>  Name  of  this  plot  in  not  more  than 

3  characters. 

2  (printer  file'  «>  Number  of  characters  per  line 

( integer  1 . 

3  (tektronix)  »>  The  baud  rate  (inleqer). 

id  -  In  less  than  30  chcrac'err.  this  is  the  box  and  id. 

mapset  -  •  0  »>  A  1  inear-i  ineci-  l  (default). 

1  =  >  A  linear-log  p  <  linear,  y  log). 

2  =*>  A  log- linear  plot  ty  1  inear.  x  log). 

3  =*>  A  log- log  plot . 

xsize  -  Tne  length  (in  inches;  of  the  x-axis  (defaults  to  6.0). 

DISSPLA  conventions  for  "inches"  are  employed,  i.e.  the 
"page"  is  assumed  to  be  9.5x11  inches. 

ysize  -  The  length  (in  inches)  of  the  y-axis  (defaults  to  4.5). 

grids  -  ■  0  »>  No  grid  (default). 

1  ■»>  Grid. 

x  -  The  user  can  also  input  the  values  of  x  from  tty  or  from 

the  input  file. 


-8- 


The  program  structure  now  loot's  as  follows: 

40  coni inue 

call  <the  subroul ine> 

if  (  nol  Ihe  first.  (  ime  ana  plotmode  >  0  )  call  <store  y> 
if  (  1(30). It. 0  )  go  to  70 
call  <muller  to  find  roots> 
write  <resu 1 t  > 
go  to  40 
70  coni  inue 

if  (  plot  desired  )  call  <plol> 
call  exit  (  abug  ) 

The  graphics  library  used  is  DISSPLA.  fill  warning  or  error 
messages,  if  any.  will  be  sent  to  the  file  named  "disout".  The 
maximum  number  of  roots  which  can  be  plotted  is  25;  however,  since 
DISSPLA  has  only  15  different  symbols  to  distinguish  different 
curves,  the  user  is  advised  not  to  plot  more  than  15  roots. 

The  plotted  output  includes: 

(1)  The  contents  of  the  input  file,  if  there  is  one. 

(2)  The  real  part  of  the  roots  vs.  the  parameter  (x) . 

(3)  The  imaginary  part  of  the  roots  vs.  the  parameter  (x) . 

The  size  of  a  page  is  assumed  to  be  8.5"xll",  regardless  of 
output  medium  selected.  If  possible.  SOLVER  will  place  items  (2) 
and  (3)  on  one  page.  The  maximum  size  of  graphs  to  be  plotted 
two-to-a-page  is  6"x4.5"  (or  4.5“x6"  for  side-by-side  plots). 


For  an  fr80  file,  the  user  must  do  either  of  the  following 
after  finding  out  the  name  of  the  f  r80  file: 

(1)  netplot  Cusc]  <fr80  file>  Cl.]  [turn.]  /tv  or 

(2)  nelout  a  <fr80  file>  b.  /tv 

(then  log-off  the  CRAY- 1  and  log-on  to  the  7608  and  view  the 
plot  file  on  a  teKtronix  by  FR80PLOT) . 

These  operations  (without  the  time  and  value)  can  be  done  on  the 
same  suffix  under  which  SOLVER  is  running. 


Example  IV. A. 1  (see  Fig.  IV. A. 1): 

A  "back-door"  way  to  plot  sin(x)  using  SOLVER 


user 

rout ine 
rout  ine 
rout ine 
rout ino/user 
rout ine/user 
rout ine/user 
rou  t ine/user 
rout ine/user 


so  1  ver 

root  solver  79.08  *** 

extracting  solv/b/1 
c  07/22/79  20:19:4!  001222 
>suh 

:  equivalence  (  ixind.l(29)  ) 

:  data  1(30)  /21/.  ixind  /0/ 

:  data  pi  /3. 1415926535/ 

:  c 
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-2.0  0.0  2.0  4-0  6.0  8.0  10.0  12.0  M.O 


PARAMETER  X 


in 


-2.0  0.0  2.0  4.0  6.0  8.0  10. 0  12.0  M.O 


PARAMETER  X 

F iqure  IV. 8. 1 


rout 

ine/user 

:  1 (30)  =  1 (30)  - 

1 

rout 

ine/user 

:  ixincf  =  ixind  + 

1 

rout 

ine/user 

:  x <  i x ind)  =  . 2  * 

0  i 

*  (  ixind-l  ) 

rout 

ine/user 

:  / 

rout 

ine/user 

>f 

rout 

ine/user 

:  equivalence  (  ixind 

.1(29)  ) 

rout 

inQ/USGT 

:  c 

rout 

ine/user 

:  f=z-sin(x( 

i  x  i 

nd)  ) 

rout 

ine/user 

>go  tty  end 

rout 

ine 

compiling  and  loading 

rout 

ine 

FT004  -  CFT  VERSION  - 

07 

/14/79  SCHEDULER 

rout 

ine 

FT001  -  COMPILE  TIME  = 

0.0130  SECONDS 

rout 

ine 

***  cray  loader  version 

- 

c 12 1  07/05/79 

rout 

ine 

exscut ing 

rout 

i  ne/user 

-p!otmode=!  $ 

rout 

ine 

roots : 

rout 

ine 

(  0.000000000000-01 

, 

0 . 000000000000-0 1 

) 

rout 

ine 

roots: 

rout 

ine 

(  5 . 877852522?8e-0 1 

, 

0.000000000000-01 

) 

rout 

ine 

roots : 

rout 

ine 

(  9.5105651 6284©-0 1 

, 

0 . 00000000000e-0 1 

) 

rout 

ine 

roots: 

rout 

ine 

(  9 . 5 10565 163 12e-0 1 

0 . 000000000000-0 1 

) 

rout 

ine 

roots : 

rout 

ine 

(  5 . 8778525235 l©-0 1 

p 

0 . 000000000000-0 1 

) 

rout 

ine 

roots : 

rout 

ine 

(  3 . 93949606966©- 1 1 

-1.89326617253e-29 

) 

rout 

ine 

roots : 

rout 

ine 

(  -5 . 87735252205e-0 1 

p 

6.31088724177e-30 

) 

rout 

ine 

roots: 

rout 

ine 

(  -9 .5 10565 16256e-0 1 

-2.24207754292e-44 

) 

rout 

ine 

roots : 

rou  t 

ine 

(  -9.51 0565 16340e-0 1 

-4.084169905270-53 

) 

rout 

ine 

root  s : 

rout 

ine 

(  -5.87785252423©-0! 

p 

-4. 21687917729 e-31 

) 

rout 

ine 

i  ools: 

rout 

ine 

(  -1  .?9599921393e-10 

p 

4.25795984001-109 

) 

rout 

ine 

roots: 

rout 

ine 

(  5. 87785252 133e-01 

p 

-1 .51273121674-123 

) 

rout 

ine 

roots: 

rout 

ine 

(  9.51056516229e-01 

p 

1 .07486017721-137 

) 

rout 

ine 

roots : 

rout 

ine 

(  9. 51056516367©-01 

p 

1.58613441594-146 

) 

rout 

ine 

roots : 

rout 

ine 

(  5.87785252496o-01 

7. 1  1282799835-161 

) 

r  out 

ine 

roo l s : 

rout 

ine 

(  2 . 69356460379e- 10 

p 

-1 .26217744835e-29 

) 

rout 

ine 

roots : 

rout 

ine 

(  -5 . S7785252060e-0 1 

p 

-1 .262177448350-29 

) 

rout 

ine 

roots: 

rout 

ine 

(  -9.510565162010-01 

p 

-1 .274473528910-57 

) 

rout 

ine 

roots: 

rout 

ine 

(  -9.510565163950-01 

p 

1.353351861820-71 

) 

rout 

ine 

roots: 
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0 . 00000000000e-0 1  ) 


routine  :  (  -5.87785252569e-01  - 

rout  ine  :  roots : 

rout  ins  :  (  -3 . 59 179842786e- 10  ,  0.00000000000e-01  ) 

rout  ins  :  all  done 

Some  users  rnay  prefer  the  following  equivalent  forms  for 
sub  and  f : 

( 1 )  sub  - 

data  1(30)  /21/.  pi  /3. 1415926535/ 

1(30)  =  1(30)  -  1 

x(21-l (30) )  =  .2  *  pi  *  (  20-1(30)  ) 

(2)  f  - 

f  =  z  -  sin  (  x(21-l(30))  ) 

The  user  may  specify  "of =3"  to  suppress  the  listing  of  roots. 


Example  IV. 8. 2  (see  Fig.  IV. ft. 2): 

ftnother  “back-door"  function  plot,  on  log-log  scales  this  time 


user 

rout ine 
rout ine 
rout ine 
rout ine/user 
rout ine/user 
rout ine/user 
rout  ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine/user 
rout ine 
rout  ine 
rout ine 
rout ine 
rout ine 
rout ine/user 
rout ine/user 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 
rout ine 


solver  sub  f  go  tty  end 
***  root  solver  79.08  *** 

extracting  solv/b/1 
c  07/22/79  20:19:41  001222 

:c  note  that  we  want  a  log  plot,  so  we  do  not 
:c  plot  points  with  negative  values  of  the 
:c  parameter. 

:c 

:  equivalence  (  ixind.l(29)  ).  (  xnext»l(28)  ) 

:  data  1(30)  /ll/»  ixind  /-&/ 

:c 

:  1 (30)  -  1 (30)  -  1 

:  plot  mode  =  ixind 

:  if  (  ixind. gt.0  )  x(ixind)  ■  xnexl 

:  ixind  =  ixind  +  1 

:  xnext  =  ixind**3 

:/ 

:  equivalence  (  xnext. 1(28)  ) 

:c 

:  f  »  z  -  xnext 

compiling  and  loading 

FT0O4  -  CFT  VERSION  -  07/14/79  SCHEDULER 

FT001  -  COMPILE  TIME  =  0.0153  SECONDS 

***  cray  loader  version  -  cl21  07/05/79 
execut  ing 

-plotmode*!  mapset=3  xsize=4.5  ysize=6.0 
-pi ot  id* " test ”  id  =  "b22"  "Stephen  au-yeung"  $ 
roots: 

(  - 1 . 25000000000e+02  .  1 .37001788954e-24  ) 

roots: 

(  -8 . 40900000000e+0 1  ,  -3.0 182941 3028e-50  ) 

roots : 

(  -2 . 70000000000e+0 1  .  -3. 1 1734638189e-67  ) 

roots: 
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rout ine 

(  -e. 000000000006+00 

7.67045853953e-93 

) 

rout ine 

roots : 

rout ine 

(  -1.00000000000e+00 

6. 16896836431-110 

) 

rout ine 

roots: 

rout ine 

(  3 . 02922587605e-28 

-1.61229026582-137 

) 

rout ine 

roots: 

rout ine 

(  1 . 0000000B000e+00 

2.17066234129-165 

) 

rout ine 

roots: 

rout ine 

(  8.00000000000e+00 

-1.75344747921-191 

) 

rout ine 

roots: 

rout ine 

(  2.70000000000e+01 

-3.39941662434-218 

) 

rout ine 

roots: 

rout ine 

(  6.400000000006+01 

5.63600517340-231 

) 

rout ine 

roots : 

rout ine 

(  1.250000000006+02 

-3.99599503859-256 

) 

rout ine 

all  done 

user 

f  i  \  es  f 

rout ine 

1375  rwe  f  105ro0x 

rout ine 

al  I  done 

user 

netplot  use  fl05ro0x  1. 

box  b22  Steve 

rout ine 

destination:  ueb 

rout  ine 

4  frames  of  output 

processed 

rout  ine 

f 105ro0x  ;  file  id  is: 

sieve  -  rxucb/rz 

rout ine 

all  done 

Example  IV.fi. 3  (see  Fig.  IV.fi. 3).  by  Y.J.  Chen: 

In  a  lei  particle-hybrid  simulation  of  the  lower — hybrid  drift 
instability  (LHDI)  in  a  uniform  magnetic  field,  the  quiet  start 
Maxwellian  loader  was  used  to  set  up  the  warm,  unmagnetized  ion 
particles.  Dispersion  relations  for  the  LHDI  and  the  multi-beam 
instabilities  were  obtained  from  SOLVER.  Initial  guesses  for  roots 
were  made  in  order  to  determine  the  approximate  frequencies  at 
which  the  desired  roots  should  be  obtained. 


user 

solver  al ib  lhdmbi(LF)f  disper  sub  mbilhd  go  imbilhd 

rout ine 

***  root  solver  79.03  *** 

rout ine 

extracting  solv/b/1 

rout ine 

c  07/22/79  20:19:41  001222 

rout ine 

compiling  and  loading 

rout ine 

FT004  -  CFT  VERSION  -  07/14/79  SCHEDULER 

rout ine 

FT001  -  COMPILE  TIME  *  0.0129  SECONDS 

rout ine 

***  cray  loader  version  -  cl21  07/05/79 

rout ine 

execut ing 

rout  ine 

routba  is  the  output  file 

rout  ine/user 

>netplot  fl05ro0x  l.  box  bll  yjehen 

rout ine 

destination:  ueb 

rout ine 

6  frames  of  output  processed 

rout ine 

f  105ro0x  :  file  id  is:  yjehen  -  rxucb/la 

rout ine/user 

>end 

rout ine 

ail  done 

The  contents  of  the  various  files  are: 
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Figure  IV. ft. 3  (a) 


(1)  disper  - 

real  v(l). lambda 
complex  t'sum(l) 

equivalence  ( v( 1 ) . p ( 1 . 33) ) .  ( f sum( 1 ) . p ( 1 , 49) ) 
equivalence  (xel.c(8)).  (xe2u. c (9) ) .  (vek.c(18)) 
equivalence  (ngr.l(2)) 

c 

f  sum(  1 )  e'-ngr*(  1 .  +><e  l+xe2u/(z-vok) ) 

ngr2-ngr/2 

do  10  i=2.ngr2+l 

f  sum( i ) =  f sum( i- 1 ) +1 . /( (z-v( i- 1 ) )*(z-v( i- 1) ) ) 

10  continue 

do  20  i =ngr2+2,ngr+l 

f  sum(  i )  =  f  surnt  i-  1 )  + 1 .  /( (z-v(  i- 1 ) )  *(z-v(  i- 1 ) ) ) 

20  coni  inue 

f  °f  sum(ngr+l ) 

(2)  mbilhd  - 

c  appears  as  the  subroutine  "vvsetup"  in  the  "solver", 
c  calculate  the  dispersion  of  the  lower  hybrid  drift  instability 
c  including  effects  of  the  multibeam  instability  due  to  the  quiet 
c  start. 

real  xx(l).vx(l).v(l).k.kmin.klast 
real  Ini.  1 l i .  1 s i 


real  lambda 

©qu 

i  val ence 

( 

XX ( 1 ) . p ( 1 .  1 )  ) 

©qu 

ival ence 

( 

vx ( 1 ) . p ( 1 . 1 7 )  ) 

©qu 

ival ence 

( 

v( 1 ) . p ( 1 . 33)  ) 

©qu 

ival ence 

( 

vt i , c( 1 )  ) 

©qu 

ival ence 

( 

vd . c ( 3 )  ) 

equ 

ival ence 

( 

upc.c(4)  ) 

©qu 

ival ence 

( 

vle.c(5)  ) 

©qu 

ival ence 

( 

ve . c ( 6 )  ) 

©qu 

ival ence 

( 

rm.c(?)  ) 

©qu 

ival ence 

( 

xel.c(8)  ) 

©qu 

ival ence 

( 

xe2w.c(3)  ) 

©qu 

ival ence 

( 

vek . c ( 1 0 )  ) 

©qu 

ival ence 

( 

b.c(ll)  ) 

©qu 

ival ence 

( 

dk . c ( 1 2 )  ) 

©qu 

ival ence 

( 

krnin.c(13)  ) 

©qu 

ival ence 

( 

Klast.c(14)  ) 

equ 

ival ence 

( 

Ini  .c( 15)  ) 

©qu 

ival ence 

( 

lti.c(lS)  ) 

equ 

ival ence 

( 

1  S  i  .  c  (  1 7  )  ) 

©qu 

ival ence 

( 

re.c(18)  ) 

©qu 

ival ence 

( 

nn, 1 ( 1 )  ) 

equ 

ival ence 

( 

ngr .1(2)  ) 

©qu 

ival ence 

( 

npts  .1(27)  ) 

c 

iff  l(29).ne.0  )  go  to  28 
1  ( 30 ) * (K 1 asl-km i n) /dk+ 1 . 5 
npts  *  1(30) 

wr  i  le (6.  1 )  nn.nqr . vt i . wpc. vt e. rm> dk . km in. k 1 ast . 
.  1  n  i ,  I  t  i .  1  s  i 
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1  format  (//."  ion:  nn-“.il5."  ngr3".il5."  vt  i  - " ,  e  14. 5. /, 

.  “el eclron : " . // , 

.  "  upc"“»el4.5. "  vle“".el5.5»/. 

.  "  mass  ratio,  rm- " . e22 . 5. /V/, 

.  “  dk-".el5.5. ”  kmin= " ,  e  1 4 . 5. *  kl ast ■ " . e 12.5. /. 

“  1  n  i  °  .  e  1 4 . 5 . "  ! t  i=".el5.5.  ’  1  s i  =  ” . e 15 . 5. ///) 

c 

c  set  op  the  multi beam 
vmax-5 . *vt i 
dv“vmax//(nn-  1 ) 
xx(  1) =0. 
do  10  i-2.nn 
w-  <  i  —  1 . 5)  *dv 
f  u-exp (- . 5*( uv/ut i ) #*2) 

10  xx ( i ) «xx (  i- 1 )  +f  u 
c 

df  *xx (nn) /ngr 

i 1 *ngr/2+l 
12*11-1 

do  12  i«l.ngr.2 
f  u- i*df 

13  if(  f v. 1 t . xx ( j+1 )  )  go  to  14 

j-j  +  1 

if(  j.gt.nn-1  )  go  to  25 
go  to  13 

14  w=du*(  j-l+(  fu-xx( j)  )/(  xx ( j+l)-xx(j)  )  ) 
vx ( i  1 )  =w 

ox (12) — w 
i 1-i 1+1 
12  12-12-1 
25  cont  inue 
c 

ngr2-ngr/2 

c  change  ute  to  the  electron  larmor  radius 

27  cont inue 

uce  i  -sqr  t  (uipc/rm) 
uc  i  i  -uce  i*rm 
re-vle*uce i 

ve-ol i *vt i#( lni+1 t i ) *uc i i 
ur i te (6. 100) uc i i .uce i ,  re. oe 

100  format(“  uc  i  i  - " .  e  13 .5.  “  ucei  =  ".el3.5. "  re*".el6.5»/'» 

.  “  ve -  * . e 1 5 . 5. //) 

28  cont inue 

1 (29)  -  1 (29)  +  1 
1 (30)  -  1 (30)  -  1 
if  (  l(30).lt.0  )  go  to  36 
c 

c  calculate  coefficients  of  the  dispersion  relation 
k-klasl-1 (30)*dK 
ur i l e ( 6 . 29 )  k 

29  format (/V. "k»a,ol2. 5) 
x ( 1 (29))  -  k 

do  30  i *  1 . ngr 
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r 


v(  i ) =VX  C i ) WK 
30  cont inue 
vek =ve*k 
bK=vte*k 
b=bk*bk 
mm=n 
c 

if  C  b.eq.0.  )  go  (.0  32 
b i0=besse i (0 . b) 
b  i  1  =besse  i  ( 1 , b) 

1 ambda=b i0*exp (-b) 
a=(  1.- lambda  ) /b 
fle-b*(  l.-bil/bi0  ) 
go  to  33 

32  lambda=1.0 
a=  1 . 0 

f le=0.0 

33  xel=a*upc 

rli=  Ini  -  lsi*(l.-fle)  -  iti*fle 
xe2u=upc*l  arnbda*r  1  i/(k>i<ucei) 
c  guess  the  roots  for  the  first  time  only 
if  (  1 (29) .gt .  1  )  go  to  36 
c 

j 

35  cont inue 

if  (  j  . g t  . miri0 (ngr2» mm)  )  return 
j 1 =ngr+l- j- 10 
y( j) =v( j l)+cmplx(0. , .001) 
j  =  j  +  l 
go  t  o  35 
c 

36  cont  inue 

do  40  j  =  1 , mm 

y  ( j )  =crnol  x  (real  (y  ( j ) ) .  abs (aimag (y ( j) ) ) ) 

40  con l  inue 

(3)  imbilhd  - 

n=4  maxit=300  h=. 00002  of =2 

1 ( 1) =16384  l (2) =16384 

c(l)=. 141421  c (3) =0 .0  c (4) = 1 . 0 

c (5) =0 . 00  c (6) =0 . 0  c  (?)  = 1600 . 0  c( 12) =0.5 

c( 13) =3.0  c ( 14) =7 . 0 

c ( 15) * . 075  c ( 16) =0 . 0  c ( 1?) =0 . 0 

y(l)-(,002. .001) 

plotrnode=l  pi ot  id= "  1  hdi "  id="box  to  1 1 "  "y  j  chen" 
xsizo=9.  ysize=6.?5 

$ 

l(l)«nn»  l(2)=ngr  for  giving  y(j)  in  input  file 
c(l)=vt»  c(3)=vd»  c(4)=upc 
c(5)=vte,  c(6)=ve.  c(7)=rm»  c(12)=dk 
c(13)=kmin»  c ( 14) =kl ast 
c(15)=lni,  c ( 16) =1 1 i »  c ( 1 7) = 1 s i 
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B. 


Plotting  a  Function  us.  One  Argument 


SOLVER  can  also  plot  a  function  vs.  x,  the  real  array 
described  in  Section  IV. A.  To  tell  SOLVER  that  the  user  wants  a 
plot  of  the  function  rather  than  the  roots  of  the  function,  the  user 
need  only  assign  the  value  4.  5.  or  6  to  plotmode  from  the  terminal 
or  from  the  input  file.  The  destination  of  the  output  when  plotmode 
equals  4.  5.  or  6  is  the  same  as  that  when  plotmode  equals  1.  » 

2.  or  3.  respectively.  All  other  variables  but  "of"  retain  the  same 
meaning  and  usage.  The  only  output  is  a  plot  in  this  case. 


Also  note  that  within  the  function  f,  z  is  equivalent  to  x(i) 
where  i  corresponds  to  the  number  of  times  the  function  has  been 
called.  The  structure  can  be  viewed  as  follows: 

40  cont  inue 

call  <the  subrout ine> 

if  (  not  the  first  time  )  call  <store  y(l)> 
if  (  1(30). It. 0  )  go  to  70 
i  =  i  +  1 
y ( 1 )  -  f  (  x(i)  ) 
go  to  40 
70  coni  inue 

call  <plot> 
call  exit  (  dbug  ) 

Example  IV.B.l  (see  Fig.  IV.B.l): 

user:  solver  f  fl/s  sub  sl/s  go  gl/d  end 

routine  :  ***  root  solver  79.08  *** 

routine  :  extracting  solv/b/1 

routine  :  c  07/22/79  20: 19:41  001222 

routine  :  compiling  and  loading 

routine  :  FT004  -  CFT  VERSION  -  07/14/79  SCHEDULER 

routine  :  FT001  -  COMPILE  TIME  =  0.0129  SECONDS 

routine  :  ***  cray  loader  version  -  cl21  07/05/79 

routine  :  executing 

routine  :  all  done 

The  contents  of  the  various  files  are: 

(1)  fl/s  - 

equivalence  (  ixind.l(29)  ) 
c 

f  =  s in  (  x ( ix ind)  ) 
c 

c  the  above  lines  for  the  function  f  are  equivalent 
c  to  the  following  alternative  line 
c  f  *  csin  (  z  ) 


REfiL(F)  IMRG(F) 

-1.5  -1.0  -0.5  0.0  0.5  1.0  1.5  -10.0  -5.0  0.0  5.0 


r 


(2)  sl/s  - 

equivalence  (  ixind. 1(29)  ) 
data  1(30)  /21/.  ixind  /0/ 

c 

1  (30)  =■  1  (30)  -  1 
ixind  ■  ixind  +  1 
x(ixind)  =  .2  *  (ixind-1  ) 

(3)  g 1/d  - 
1  (30) »6  1 

plotmode=4  grids=l 
of  =3 
* 


C.  Plotting  Roots  (a  Function)  vs.  Two  Parameters  (Arguments) 


A  pair  of  contour  plots  will  be  produced  in  this  case.  The 
second  parameter  (arqument)  is  to  be  stored  in  w.  a  real  array  of 
size  1024.  The  user  should  store  values  into  x  and  w  in  the 
subroutine.  In  addition  to  x  and  w«  two  more  variables  have  to 
be  set  at  execution  time:  they  are: 

nx  -  Humber  of  elements  in  x  (defaults  to  1024). 


nw  -  Number  of  elements  in  w  (defaults  to  25). 


If  the  product  of  nx  and  nw  is  greater  than  1024x25=25800. 
nx  and  nw  will  be  set  bacK  to  their  default  values.  It  is  important 
that  w  be  the  more  rapidly  varying  parameter  of  the  two  (x  is  the 
slowly  varying  parameter).  This. point  will  become  clear  upon 
examination  of  the  examples  which  follow. 


Example  IV.C.l  (see  Fig.  IV.C.l): 

In  this  example,  the  fastest  growing  "tilting  mode"  root  of 
a  weaK  ion  ring  -  plasma  system  is  identified  and  contours  of  its  real 
frequency  and  growth  rate  normalized  to  the  ion  cyclotron  frequency 
are  plotted  as  functions  of  two  parameters.  These  are  ths  ring 
strength  "qr"  and  the  background  plasma  density  "zni".  .garithmic 
scales  are  used,  and  (in  the  stable  region  at  the  lower  right)  the 
frequency  is  set  to  zero  when  the  growth  rate  is  essentially  zero, 
since  there  is  then  no  "fastest  growing  mode".  Note  that  the  roots 
are  reordered  since  the  contour  plotter  can  only  plot  one  root,  and 
we  have  not  arrcnged  to  "follow"  individual  roots  as  parameters 
change.  The  constant  “cay"  depends  upon  the  differencing  scheme 
being  used.  (A.  Friedman,  unpublished) 

user:  solver  f  f6/s  sub  s6/s  go  g6/d 

routine  :  *#*  root  solver  79.08  *** 

routine  :  extracting  solv/b/1 
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ROO 


IMflG (ROOT) 


PARAMETER  X 

Figure  IV.C.l  (b) 


rou lino 
roul ino 
rout ino 
rout ino 
rotu ine 
roul ino 
rout  inc?/user 
rout ine 
rout ine 
rout ine 
rout ine 
roul  ineAjser 
rout ine 


c  07/28/79  11:21:15  00122? 
compiling  and  loading 

FT004  -  CF-T  VERSION  -  07/14/79  SCHEDULER 

FT0LI1  -  COMPILE  TIME  =■  0.0245  SECONDS 

**#  cray  loader  version  -  cl?l  07/05/79 
oxocut ing 

>netplot  f 105ro0x  1. 
box  number :  b22 

destination:  ucb 

6  frames  of  output  processed 
f 105ro0x  ;  file  id  is:  fl05ro0x  -  rxucb/4K 
>end 

al 1  done 


The  contents  of  the  various  files  are: 


(1)  f 6/s  - 

equivalence  C  qr.l(26)  ).  (  cay. 1 (27)  ).  (  zni.l(28)  ) 
c 

f  =  z  *  z  *  (  l.+z  )#*2  -  qr  *  (  z*z/3 .  +2 .  *cay/zn  i  ) 

(2)  s6/s  - 

equivalence  (  qr.l(26)  ) ..  (  cay.  1(27)  ),  (  zni.l(28)  ) 

equivalence  (  ipoint.l(29)  ) 

data  cay/. 25/.  ipoint/0/.  1(30)/-1000/ 

c 

if  (  1 (30) .eq.-1000  )  1(30)  =  rx  *  nu 

1 (30)  =  1 (30)  -  1 
i point  *  i  point  +  1 
c 

ixind  =  (  ipoint-1  )  /  nx  +  1 

znilog  =  (  ixind-1)  /  10. 
zni  =  1 0 . **zn i 1 og 
x (  ix  ind)  =  zn i 1 og 
c 

iuind  =  mod  (  ipoint-1.  nu  )  +  1 
qrlog  =  (  iuind-1  )  /  10. 
qr  •  .36  *  10.**qrlog 
u( iwind)  =  qrlog 
c 

c  sort  the  roots  so  that  the  fastest  growing  mode  appears  first 
n  =  4 

c(l)  =  y(l) 

do  2  iroot  =2.  4 

if  (  a  irnag  (y  (  iroo  t ) )  .  1  e .  a  imag  (c  ( 1 ) )  )  go  to  2 
c(  1 )  =  y(  iroot) 
y( iroot)  =  y ( 1) 
y ( 1 )  =  c(l) 

2  cont inue 
c 

c  set  root  to  zero  if  growth  rate  is  below  a  reasonable  threshold 
if  (  a imag(y ( 1 ) ) . 1 t . 0. 0 1  )  y(l)  =0. 

(3)  g6/d  - 

nx“31  nu=31 

pi ot mode*  1  mapset=4 
xsize-6.  ysize=6. 
of  =3 
$ 
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Example  IV. C. 2  (see  Fig.  IV. C. 2): 


user 

rout ine 
rout ine 
rout ine 
rout  ine 
rout ine 
rout ine 
rotuine 
rout ine 
rout ine/user 
rout ine 
rout  ine 
rout ine 
rout ine 
rout ine/user 
rout  ine 


solver  f  f5/s  sub  s5/s  go  g5/d 
***  root  solver  79.08  *** 

extracting  solv/b/1 
c  07/28/79  11:21:15  001222 
compiling  and  loading 

FT004  -  CFT  VERSION  -  07/14/79  SCHEDULER 

FT001  -  COMPILE  TIME  »  0.0245  SECONDS 

***  cray  loader  version  -  cl21  07/05/79 
execut ing 

>netplot  fl05ro0x  1. 
box  number:  b22 

destination:  ucb 

6  frames  of  output  processed 
fl05ro0x  ;  file  id  is:  f  105ro0x  -  rxucb/\; 
>end 

all  done 


The  contents  of  the  various  files  are: 


(1) 

c 


f  5/s  - 

equivalence  (  ixind. 1(29)  ).  (  iuind.  1(28)  ) 
f  -  0 

rsq  =  (  abs(x( ixind) )-.5  )**2  +  . 25*u( iuind) **2 
if  (  rsq. ge. 0.25  )  return 
f  =  cmplx  (  sqrl ( ,25-rsq) .  sqrt(rsq)  ) 
if  (  x( ixind) .gt .0.  )  f  =  -f 


(2)  s5/s  - 

equivalence  (  ixind. 1(29)  ).  (  iuind. 1(28)  ) 
equivalence  (  denx.l(27)  ).  (  denu. 1(26)  ) 
data  1 (30)/-1000/.  ixind/0/.  iuind/25600/ 

c 

if  (  1 (30) . ne . -  1000  )  go  to  10 
1  (30)  =  nx  *  nu 
denx  =  nx  /  2 
denu  =  nu  /  2 
c 

10  cont inue 

1  (30)  =  1 (30)  -  1 
if  (  iuind.lt.nu  )  go  to  20 
ixind  =  ixind  +  1 

x(ixind)  =  1.2  *  (  ixind-l-denx  )  /  tienx 
iuind  =  0 
20  cont inue 

iuind  =  iuind  +  1 

u(iuind)  =  1.2*  C  iuind-l-denu  )  /  denu 


(3)  g5/d  - 
pi otmode=4 
mapsel =4 
nx«81  nu*81 

xsize=7.  ysize=7. 
of  =3 
$ 


-26- 


-0.9  -0.6  -0.3  0.0  0.3  0.6 

ARGUMENT  X 


Ref  prences 


1.  S.D.  Conte  and  C.  d©  Boor,  "Elementary  Numerical  Analysis: 

An  Algorithm  Approach."  Ed.  2.  pp.  7 4-83.  McGraw-Hill. 

New  York  C 1972) . 

2.  M.J.  Geruer,  C.K.  Birdsall.  A.B.  Lanqdon  and  D.  Fuss.  "Normal 

Modes  of  A  Loss  Cone  Plasma  Slab  with  Steep  Density  Gradient 
Memorandum  No.  ERL-M602.  Electronics  Research  Laboratory, 
College  of  Engineering,  UC  Berkeley.  <24  September  1976). 

See  Appendix  I  for  a  description  and  listing  of  ROOTS  which 
includes  MRAF  and  FOLLOW,  which  traces  roots,  including 
branch ing. 

3.  M.J.  Gercer.  "Roots.  A  Dispersion  Equation  Soluor,"  Memorandum 

No.  ERL-M77-27.  Electronics  Research  Laboratory,  College  of 
Engineering,  UC  Berkeley  (31  October  1976). 

This  report  complements  Ret .  2  and  contains  instructions  to 
the  user  of  ROOTS  (MRAF,  FOLLOW,  etc.). 

4.  “DISSPLA  User's  Manual."  Version  8.0.  ISSCO.  San  Dieqc  (1978). 

5.  "Third  Quarter  Progress  Report  on  Plasma  Theory  and  Simulation. 

ERDA  Contract  FT-  c-S-03-0034.  Project  Agreement  No.  128. 
Electronics  Research  Laboratory.  Colleqe  of  Engineering. 

UC  Berkeley  (1  October  1977) . 


-29- 


APPENDIX  I. 


The  Structure  of  The  Library  Solver 


As  mentioned  in  Section  I,  SOLVER  is  a  1  ibrary  that  contains 
all  files  needed  by  the  SOLVER  program. 


Entering  "solver  /  t  v"  from  the  terminal  causes 
"solver'x",  the  first  entry  in  SOLVER.,  to  be  executed.  A  dropfile 
named  +solver<s>.  where  <s>  is  the  current  suffix  (channel)  that 
SOLVER  is  running  under,  is  created.  The  next  step  is  extracting 
the  binary  library  "solv'b'1"  from  "solver"  (if  "solv'b'1" 
is  not  yet  in  the  user's  active  file  area).  The  prompt  ">“  will 
then  appear  and  the  user  proceeds  to  enter  commands. 


The  rest  of  the  files  in  the  SOLVER  1 ibrary  are:  the  latest 
version  of  this  documentation  and  the  source  files.  Some  of  these 
files  are  part  of  the  SOLVER  controller  and  the  rest  are  used  in 
the  generated  program.  A  list  of  routines  and  correspond  ing 
source  files  follows: 


(1)  The  SOLVER  controller  program  - 
file  name  routine  name 


sol ver's 

main,  (the  mi 

cmds's 

iscrnd 

digit's 

digit 

f  ex ist's 

f  ex ist 

f  unct's 

get  f ct 

inHe's 

i  n  i  t  e 

letter's 

1 et  ter 

prog's 

putprog 

readf 's 

readfct 

remove's 

remove 

savef 's 

savef 

scncmd's 

scancmd 

sf msge's 

sh i f  tmsg 

symove's 

symove 

»  SOLVER  generated  program  - 

file  name 

rout ine  name 

dr iver's 

driver  (the  i 

genei 

curves's 

curves 

dest's 

dest  in 

map l yp's 

maptype 

root's 

rnul ler 

ncyc's 

ncyc 

plot's 

rootpl ot 

se l f 's 

setf  ile 

step's 

stepsz 

wse  t 's 

wse  tup 
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To  modify  SOLVER,  the  user  need  only  recompile  those  source 
files  that  are  modified  and  use  BUILD  to  update  "sol v/b/1 " .  If 
all  the  modified  files  are  in  item  (2)  above.  the  user  need  only 
update  the  SOLVER  1 ibrary  “solv/b/l";  otherwise,  the  user  must 
generate  a  new  "sol ver/x'  by  using  LDR  as  follows: 

1  dr  i =sol ver/b.x  -sol ver/x. 1 ib=sol v/b/1 : f or  1 1 ib 

"solver/b"  is  the  binary  obtained  by  compiling  main,  (in  file 
“solver/s")  using  CFT.  If  “solver/s"  is  not  modified,  ''solver/b" 
can  be  generated  by: 

build  ol .  soi v/b/1 (LF)g.  sol ver/b(LF)xg.  main. (LF lend 
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APPENDIX  I.  The  Structure  of  The  Library  Solver 

As  mentioned  in  Section  I.  SOLVER  is  a  library  that  contains 
all  files  needed  by  the  SOLVER  program. 

Entering  "solver  V  t  v“  f<-om  the  terminal  causes 
"sol ver/x " .  the  first  entry  in  SOLVER*  to  be  executed.  A  dropfile 
named  +solver<s>.  where  <s>  is  the  current  suffix  (channel)  that 
SOLVER  is  running  under,  is  created.  The  next  step  is  extracting 
the  binary  1 ibrary  “sol vVbVi  “  from  “solver”  (if  "solv/b/1" 
is  not  yet  in  the  user's  active  file  area).  The  prompt  ">"  will 
then  appear  and  the  user  proceeds  to  enter  commands. 


The  rest  of  the  files  in  the  SOLVER  1  ibrary  are:  the  latest 
version  of  this  documentation  and  the  source  files.  Some  of  these 
files  are  part  of  the  SOLVER  controller  and  the  rest  are  used  in 
the  generated  program.  A  list  of  routines  and  correspond ing 
source  files  follows: 


(1)  The  SOLVER  controller  program  - 
file  name  routine  name 


sol verVs 

main.  ( t 1 

cmdsVs 

iscmd 

a ig i tVs 

digit 

f  ex  is t Vs 

f  ex ist 

functvs 

get  f  ct 

ini teVs 

in  i  te 

let  ter Vs 

1 et  ter 

progVs 

putprog 

readf  Vs 

readf ct 

remove/s 

remove 

save f Vs 

savef 

scncmdVs 

scancmd 

sf msgeVs 

sh i f  tmsg 

symoveVs 

symove 

i  SOLVER  generated  progn 

file  name 

rout  ine  i 

dr iverVs 

driver  ( 

curvesVs 

curves 

dest Vs 

dest  in 

maptypVs 

maptype 

root  Vs 

mu  1 1  er 

ncycVs 

ncyc 

pi Ol Vs 

rootpl ot 

se  t  f  vs 

set  f  i  1  e 

stepvs 

stepsz 

wset  Vs 

wset  up 

generated  program) 
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To  modify  SOLVER-  the  user  need  only  recompile  those  source 
files  that  ore  modified  and  use  BUILD  fo  update  “sol  V'b/1  “ .  If 
all  fhe  modified  files  are  in  item  (2)  above,  the  user  need  only 
update  the  SOLVER  library  “sol v/b/1 “ ;  otherwise,  the  user  must 
generate  a  new  “solver/x"  by  using  LDR  as  follows: 

1  dr  i =sol ver/b.x-sol ver/x. 1 ib=sol v/b/1 : f or 1 1 ib 

"solver/'b"  is  the  binary  obtained  by  compiling  main,  (in  file 
“solver/s")  using  OFT.  If  “solver/s"  is  not  modified,  "solver/b" 
can  be  generated  by: 

build  ol .  sol v/b/1 (LF) g .  sol  ver  Zb  (LF  ■>  xg .  main.(LF)end 
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APPENDIX  II.  Listings  of  The  Source  Files 


The  following  is  a  listing  of  all  source  programs  used  by 
SOLVER;  files  are  current  as  of  August  31.  1979,  and  are  listed 
in  the  same  order  as  in  the  table  above: 
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1  c 

2  c 

3 

4 

5 

6 

?  c 
8 

9  c 


main  program  of  solver 


i  mp  1  i  c  i  l 
common 
da  l  a 
da  l  a 


integer (a-z) 
/dabk/  Keep, 
Keep  /0/ 
termin  /0/ 


ini tprm 


integer  msg (10),  symb(4),  type (4),  len 


10 

integer  1  ibuf (20),  nl  ib,  buf (22) 

1 1 

data  buf  /"  i=#b,x",  "=#x.lib=",  "solv/b/l ' 

12 

":tv301ib",  fort  lib"/ 

13 

data  nl ib  /32/ 

14 

equivalence  (  libuf(l),  buf (3)  ) 

15 

c 

16 

c 

initial izat ion 

1? 

cal  1  dropf i 1 e  (  0  ) 

18 

call  msgflag  (  msgset,  2  ) 

19 

call  msgl ink  (59,  1  ) 

20 

call  mprompt  (  “>",  1  ) 

21 

call  rnsgtor  (  "***  root  solver  79.08  ***" 

,  32 

22 

c 

23 

c 

extract  binary  modules  from  the  library 

24 

if  (  f ex ist ( 1 ibuf ( 1 ) ) . eq . 0  )  go  to  10 

25 

call  rnsgtor  (  "extracting  solv/b/l",  20  ) 

26 

call  inite  (  “lib",  60000000,  0.,  ires  ) 

2? 

call  msgtoe  (  "solver",  6,  1,  0,  ires  ) 

28 

call  rnsgfre  (  msg,  len,  26,  ires  1 

29 

call  rnsgtor  (  msg,  len  ) 

30 

call  msgtoe  (  "x  solv/b/l",  10,  1,  0,  ires  ) 

31 

call  rnsgfre  (  msg,  len,  32,  ires  ) 

32 

if  (  len.gt.5  )  call  rnsgtor  (  msg,  len  ) 

33 

call  bypass  (  1,  1  ) 

34 

call  msgtoe  (  "end",  3,  1,  0,  ires  ) 

35 

call  rnsgfre  (  msg,  len,  16,  ires  ) 

36 

c 

3? 

10 

cont  inue 

33 

call  rnsgtor  (  "  ",  1  ) 

39 

c 

scan  commands 

40 

call  scancmd  (  "#4".  1  ibuf ,  nl ib,  termin 

) 

41 

go  to  (  20,  50,  40.  12,  40  ) ,  termin+1 

42 

c 

compile  root  solver 

43 

12 

cont  inue 

44 

termin  =  4 

45 

20 

cont inue 

46 

call  rnsgtor  (  "compiling  and  loading",  21 

) 

43 

call  inite  (  "eft",  6000000000,  0.,  ires 

) 

43 

call  destroy  (  "#bn  ) 

49 

call  bypass  (  1.  1  ) 

50 

call  msgtoe  (  " i =#i ,b=#b, 1 =0 ",  13,  0,  0, 

ires 

51 

call  rnsgfre  (  msg,  len,  80,  ires  ) 

52 

c 

load  root  solver 

53 

30 

cont  inue 

54 

call  inite  (  "Idr",  6000003000,  0.,  ires 

) 

" idisspla", 
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55 

56 
5  7 

58 

59  c 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

7 1 

72 

73 

74 

75 

76  c 

77  c 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87  c 

88 


call  destroy  (  "#x"  ) 
call  bypass  (  1,  1  ) 

call  rnsgtoe  (  buf.  nl ib+16.  0.  0,  ires  ) 
call  msgfre  (  msg.  len.  80.  ires  ) 
execute  root  solver 

call  msgtor  C  "executing".  10  ) 

40  coni  inue 

call  inite  (  11  #x 11 .  6000000000.  0..  ires  ) 
if  (  ires.ne.0  )'  go  to  45 
call  bypass  (  1.  1  ) 
nc  ■  - 1 

if  (  termin.eq.4  )  nc  =  8 
call  msgtoe  (  initprm.  nc,  0,  0,  ires  ) 
call  msgfre  (  msg,  len,  80.  ires  ) 
call  getsymb  (  symb.  type,  nc.  msg,  len.  4  ) 
if  (  syr.ib  (3)  .  ne .  “al  1  “  .or.  symb  (4)  .  ne  .  "done  “  ) 
call  msgtO'  (  msg,  len  ) 
go  to  10 
45  co.  ’.nue 

call  msgtor  (  "program  not  yet  generated11,  30  ) 
go  to  10 

destroy  temporary  files 
50  cont  inue 

call  destroy  (  "#i"  ) 

call  destroy  (  "#b"  ) 

call  destroy  (  "#x“  ) 

call  destroy  (  "i#b"  ) 

call  destroy  (  "Soul"  ) 
call  destroy  (  "disout"  ) 

if  (  keep.eq.0  )  call  destroy  (  "solv/b/1 "  ) 
call  exit 

end 
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logical  function  digit  C  ch  ) 
integer  ch 


1 
2 

3  c 

4  c  digit  returns  true  if  ch  is  a  digit 

5  c 

6  integer  zero-  nine 

7  data  zero  /  60b  A  nine  /  7 lb 

8  c 

9  digit”. true . 

10  if  (  ch.ge.zero  .and.  ch.ie.nine  ) 

11  digit  =  .false. 

12  return 

13  end 


/ 


return 
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55 

56 

57 

58 

59  c 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76  c 

77  c 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87  c 

88 


call  destroy  (  "#x"  ) 
call  bypass  (  1,  1  ) 

call  msgtoe  (  buf.  nl ib+16.  0.  0.  ires  ) 
call  msgfre  (  rnsg.  len.  80.  ires  ) 
execute  root  solver 

call  msgtor  (  "executing".  10  ) 

40  cont  inue 

call  inite  (  "#x".  6000000000.  0.,  ires  ) 
if  (  ires.ne.0  )'  go  to  45 
call  bypass  (  1,  1  ) 
nc  =  - 1 

if  (  termin.eq.4  )  nc  =  8 
call  msgtoe  (  initprm,  nc.  0.  0,  ires  ) 
call  msgfre  (  msg.  len.  80.  ires  ) 
call  getsymb  C  symb.  type.  nc.  msg.  len.  4  ) 
if  (  syr.ib  (3)  .ne .  "al  1"  .or.  symb  (4)  .  ne .  “done " 
call  msgtor  (  msg.  len  ) 
go  to  10 
45  continue 

call  msgtor  (  "program  not  yet  generated" 
go  to  10 

destroy  temporary  files 
50  cont inue 

call  destroy  (  "#i“  ) 

call  destroy  (  "#b"  ) 

call  destroy  (  “#x“  ) 

call  destroy  (  "l#b"  ) 

call  destroy  (  “$out"  ) 
call  destroy  (  “disout"  ) 

if  (  keep.eq.0  )  call  destroy  (  “sol v/b/'l  * 
coll  exit 


end 


30  ) 
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1  logical  fund  ion  digit  (  ch  ) 

2  integer  ch 

3  c 

4  c  digit  returns  true  if  ch  is  a  digit 

5  c 

6  integer  zero,  nine 

7  data  zero  /  60b  /,  nine  /  ?lb  / 

8  c 

9  digit”. true . 

10  if  (  ch.ge.zero  .and.  ch.le.nine  )  return 

11  d ig i t  =  . f al se . 

12  return 

13  end 


U)co-Ncn(_n.&.ojr\; 


i 


10 
1 1 
12 


integer  function  fexist  (  name  ) 
implicit  integer(a-z) 
c 

c  check  if  name  is  an  ex  1st ant  file 
c  fexist  returns  0  =>  f.le  exists 
c  else  =>  file  not  exists 

c 

call  f roe ioc  (  ioc  ) 

fexist  =  izopen  C  ioc*  name-  len-  acs  ) 
if  (  fexist. eq.0  )  call  izclose  (  ioc. 
return 
end 


acs  ) 
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—  OkDOD-siiTiLnx^GKjro  —  OkDCD^icnaix^CKirj  —  ouDco-Njcncn^uiro 


1  subroutine  getfct  (  buffer,  nl  ) 

integer  buf f er (9. 1 ) .  nl 
c 

c  gets  the  function  or  the  subroutine  and  stores  it  in  buffer, 
c  terminates  uhen  the  first  non-blanK  character  is  f  folloued  by 
c  a  blank  or  a  equal  sign,  or  is  a  slash 
c 

logical  digit 

integer  1 ine (80) 

integer  blank,  equ.  f.  slash 

data  blank  /  40b  / 

data  equ  /  75b  /,  f  /  146b  / 

data  slash  /  57b  / 

c 

call  mprompt  (  1  ) 

nl  =0 
20  cont inue 

nl  =  nl  +  1 

call  msgfrr  (  bufferfl.nl).  len.  72  ) 
call  zrjcharz  (  line,  buffer(l.nl).  len  ) 
i  =  0 

25  cont inue 

i  =  i  +  1 

if  (  1 ine(  i) ,eq. slash  )  go  to  40 

if  (  1 ine C i ) . eq . bl ank  .or.  d  ig  i  t C 1  ine C  i ) )  )  go  to  25 
if  C  lineCO.ne.f  )  go  to  20 

if  C  1 ineC i+1) .ne. blank  .and.  1 ine C i+1 ) . ne . equ  )  go  to  20 
c 

30  continue 

call  mprompt  C  1  ) 

return 
c 

40  cont inue 

nl  =  nl  -  1 
go  to  30 
c 

end 
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1 

subroutine  inite  (  name,  time. 

2 

i  mp  1  i  c  i  t 

integer (a-z) 

3 

integer 

name.  time,  result 

4 

real 

pr  io 

5 

c 

6 

c 

initiates  the  controllee  -  name 

7 

c 

a 

parameter 

(  nbeta=8  ) 

9 

integer 

beta(nbela) 

10 

data 

beta(l)  /20030b/ 

11 

data 

bet a (6)  /0/ 

12 

data 

beta (7 )  /0/ 

13 

data 

beta(8)  /0/ 

14 

c 

15 

real 

pr ior  i  ty 

16 

equivalence  C  priority.  beta(5) 

1 7 

c 

18 

beta(3)  = 

name 

19 

beta(4)  = 

i  imQ 

20 

priority 

=  pr  io 

21 

resul t  = 

izcray  C  beta,  nbeta  ) 

22 

return 

23 

end 
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uj  r\j  —  oi£>CD-g<T>cn.t>uif\j 


1 


c 

c 

c 


c 


logical  function  letter  (  ch  ) 
integer  ch 

letter  returns  true  if  ch  is  a  letter 
integer  a.  z 

data  a  -  141b  / .  t.  /  172b  ✓ 

letter  =  . true . 

if  (  ch.ge.a  .and.  ch.le.z  )  return 

1 et  ter  •  .false. 

return 

end 
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2 

3 

4 

5 

6 
7 


subroutine  putprog  C  ious,  ilext.  buffer,  nl .  w  ) 
implicit  integer (a-z) 

integer  ious.  ilext.  bufferO.l).  nl.  w 

writes  the  function  f  or  the  subroutine  vvsetup  onto  the  disk 
file  -  it ext 


3 

real 

h 

. 

9 

common 

7miscbk7  n 

10 

data 

n  717 

1 1 

c 

12 

integer 

fun(9),  sub(9) 

.  f ors (9. 2) 

13 

equivalence  ( 

f un( 1 ) . f ors ( 1 . 1 )  ).  (  sub(l).fors 

14 

data 

f  un 

7"  CO 

".“mplex  fu". "net  ion  " 

15 

data 

sub 

/“  su 

" .  "broul  ine " .  11  wsetu" 

16 

c 

1? 

if  (  w.gt 

.  1  ) 

go  to  10 

18 

call  freeus  ( 

ious  ) 

19 

cal!  create  ( 

ious.  itext,  2.  -1  ) 

20 

c 

21 

c 

write  program  on 

disk 

22 

wr i te ( ious. 101) 

23 

10 

cont inue 

24 

wr  i  te ( ious. 102)  (forsti 

,  w) . i  =  1 . 9) .  n 

25 

wr i te ( ious. 103)  ( Cbuf f er ( j .  i ) .  j  =  l,9),  i  =  l.nl) 

26 

wr i te( ious. 104) 

27 

return 

28 

c 

29 

101 

formal  ( 

II 

cal  1  dr i ver "  7 

30 

. 

H 

end"  ) 

31 

102 

format  ( 

9a8  7 

32 

" 

complex 

z"  7 

33 

II 

integer 

of.  id  C 10) •  plotmode. 

34 

II 

common 

7miscbk7  n.  maxit.  h. 

35 

" 

data 

n  /".  i4.  "/"  7 

36 

H 

real 

x( 1024),  w C 1024) "  7 

37 

II 

common 

/plolbk7  id.  mapset. 

38 

M 

. 

plot  id,  x,  w,  xsize.  y 

39 

N 

. 

grids,  nx.  nw"  / 

40 

II 

compl ex 

c<30) .  p(  1024.  100)  .  l) ( 

41 

II 

integer 

1 (30) 0  / 

42 

n 

common 

7parmbk/  1.  c.  p,  y" 

43 

103 

format  ( 

9a3  ) 

44 

104 

format  ( 

II 

return “ 

/ 

45 

. 

II 

end"  ) 

46 

c 

47 

end 

"f  (z)".5*"  V 


epl.  ep2.  of  *  7 


pi  ot mode. "  / 
size,  ini en.  “  / 
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1  subroutine  readfct  (  ious.  name,  fund,  n!  . 

2  implicit  integer(a-z) 

3  integer  ious.  name.  funct(9.1).  nl .  ierr 

4  data  f  irst  /l/ 

5  c 

6  c  reads  rune' ion  or  subroutine  from  a  file 

7  c 


8  if  (  f irst .eq.  1  ) 


9 

f irst  =  0 

10 

ierr  =  0 

1 1 

call  open  (  iou 

12 

if  (  ierr . 1 e  .  0 

13 

nl  =  0 

14 

10  coni inue 

15 

nl  =  nl  +  1 

16 

cal  1  rdl ine 

17 

if  (  1  on . gt . 

18 

nl  =  nl  -  1 

19 

call  close  ( 

20 

return 

2 1  c 

22 

end 

Call  freous  C  ious  ) 


name.  0.  ierr  ) 
return 


ious.  functn.nl).  len  ) 
)  go  to  10 

DUS  ) 


ierr  ) 
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1  subroutine  remove  (  name.  sir.  len  ) 

2  implicit  inleger(a-z) 

3  integer  name,  str(l).  len 

4  c 


j  c  remove  name  t rom  str.  i  en 

6  c 

7  parameter  (  cpu=8  ) 

8  parameter  (  blank=40b  ) 

9  integer  buf (25) 


10 

c 

1 1 

call  zmovechr  (  buf.  0. 

12 

c 

13 

c 

compute  the  length  of  name 

14 

c 

15 

n  *  zscanrne  (  name.  0. 

16 

c 

17 

c 

removG  name 

18 

c 

19 

1 oc  *  zskeybyt  (  buf 

20 

if  (  loc.eq.0  )  go 

21 

if  (  loc.eq.len  )  g 

22 

call  zmovechr  (  str. 

23 

len  =  len  -  n  -  1 

24 

call  zmovechr  (  str. 

25 

return 

26 

c 

27 

30 

cont  inue 

28 

n  r,  +  ) 

29 

len  =  1 en  -  n 

30 

call  zmovechr  (  sir. 

31 

return 

32 

c 

33 

c 

not 

f  ound 

34 

c 

35 

40 

coni inue 

36 

ur  i  te (59. 10 1 )  name 

37 

return 

38 

c 

39 

101 

format  (  "not  found:  “ , 

40 

c 

41 

end 

is  the  number  of  character  in  str 


str,  0,  len  ) 


cpu,  blank  )  +  cpu 


0,  len.  name,  n  ) 
o  30 
to  40 

0,  buf.  0.  loc-1  ) 

loc-1.  buf.  n+loc.  len-loc+1  ) 


0.  buf.  n.  len  ) 


a9  ) 
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1 

subroutine  savof  (  name.  fund,  nl  ) 

2 

implicit  integor(a-z) 

3 

integer  name,  fund  (9,1).  nl 

4 

c 

5 

c 

save  fund  to  the  file  "name" 

6 

c 

7 

if  (  f ex  is t (name) . ne . 0  )  go  to  5 

8 

answer  =  “n“ 

9 

call  mprompt  (  “destroy  existence 

10 

call  rnsgfrr  (  answer,  len,  1  ) 

1 1 

call  mprompt  (  1  ) 

12 

if  (  answer .ne. "y"  )  return 

13 

5 

cont inue 

14 

call  freeus  (  ious  ) 

15 

call  create  (  ious,  name.  2,  -1,  ierr 

16 

if  (  ierr.ne.0  )  go  to  20 

17 

wr  i  le ( ious. 10 1 )  ( ( f unct ( j . i ) ,  j-1,9) 

18 

call  closo  (  ious  ) 

19 

return 

20 

c 

21 

20 

cont inue 

22 

call  msglor  (  "error  creating  file 

23 

return 

24 

c 

25 

101 

format  (  9a8  ) 

26 

c 

27 

end 

24  ) 
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1 

subroutine  scancmd  (  itext.  1 ibuf .  nl ib.  termin  ) 

2 

impl ic i l 

integer (a-z) 

3 

integer 

itext.  libuf(l).  nl ih.  termin 

4 

c 

5 

c 

this  subroutine  scans  the  input  line  and  executes  the 

6 

c 

given  commands 

7 

c 

8 

parameter 

(  ncmd*10  ) 

9 

integer 

command ( ncmd ) 

IQ 

data 

command  / 

1 1 

"sub  “ . 

12 

"Kb  in". 

13 

Bdl ib”. 

14 

”xeq". 

15 

"saves “ . 

16 

“savef “ . 

17 

"al ib". 

18 

"end". 

19 

“go". 

20 

, 

"f "  / 

21 

common 

/miscbKx'  n 

22 

c 

23 

common 

/dabk/  Keep,  initprm 

24 

c 

25 

parameter 

(  maxsym=80  ) 

26 

integer 

symb (maxsym) .  t ype (maxsym) .  buffer(10) 

27 

c 

28 

integer 

i .  numsym 

29 

data 

i  /  0  /.  numsym  /  0  / 

30 

c 

31 

integer 

f unct (9. 300) .  nl 

32 

data 

funct  /  ■  f«".  “0.".  7*‘  *  / 

33 

data 

nl  /  1  / 

34 

integer 

prog (9 . 300) .  npl 

35 

data 

prog  /  "c".  8*"  "  / 

36 

data 

npl  /  1  / 

37 

c 

38 

integer 

buf  (20) 

39 

c 

40 

if  (  i.lt 

.numsym  )  go  to  20 

41 

c 

42 

c 

main  loop  to 

interpret  commands 

43 

c 

44 

10  cont inue 

45 

i  =  0 

46 

call  msgfrr  (  buffer,  len,  maxsym  ) 

47 

call  getnumb  (  symb,  type,  numsym,  buffer,  len.  maxsym  ) 

48 

typo(numsym+l )  »  7 

49 

20  cont inue 

50 

i  3 

i  +  1 

51 

go 

to  (  100*  400,  20*  420*  420.  420.  420*  10  )* 

type( i )+l 

52 

c 

53 

100  cont inue 

54 

go 

to  (  120,  300.  280.  260.  240.  230.  220.  200. 

180. 
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160.  140  ) ,  iscmdfsymbf i) .command. ncmd)+l 


55 

1 

160.  140  1.  iscmdfsymbf i) .command. ncmd)+l 

56 

c 

5  7 

c 

control  lee?  assumed 

58 

c 

59 

120 

cont inue 

60 

call  inite  C  symb(i),  60000000000,  0.,  ires  ) 

61 

if  (  ires.ne.0  )  go  to  130 

62 

cal  1  bypass  (  *1 ,  1  ) 

63 

ns  =  - 1 

64 

if  (  lype( i+1) .oq.7  1  go  to  125 

65 

call  shiftmsg  (  buffer,  len  ) 

66 

ns  =  len 

67 

125 

cont inue 

68 

call  msgtoe  (  buffer,  ns.  0,  0,  ires  ) 

69 

call  msgfre  (  buffer,  len.  maxsym.  ires  ) 

70 

call  getsymb  C  symb.  type,  numsym,  buffer,  len.  4  ) 

71 

if  (  syrnb(3)  .eq.  "al  1  *  .and.  sym.b(4)  .eq.  "done"  )  go 

72 

call  rnsgtor  (  buffer,  len  ) 

73 

go  to  10 

74 

130 

cont inue 

75 

call  rnsgtor  (  "invalid  command  d isregarded " .  30  ) 

76 

go  to  20 

77 

c 

78 

c 

f  command 

79 

c 

80 

140 

cont inue 

81 

if  f  type ( i+1 ) . ne . 0  )  go  to  155 

82 

call  roadfct  (  ious,  symb(i+l).  funct.  nl.  ires 

83 

if  (  ires. It. 0  )  go  to  155 

84 

i  s  i  +  1 

85 

go  to  20 

86 

155 

cont inue 

87 

call  getfct  (  funct,  nl  ) 

88 

go  to  20 

89 

c 

30 

c 

go 

command 

91 

c 

92 

160 

cont inue 

93 

term  in  =  0 

94 

call  putprog  C  iou.  itext.  funct,  nl.  1  ) 

95 

if  <  icall.eq.l  )  call  putprog  (  iou.  itext.  prog. 

96 

cal  1  close  (  iou  ) 

97 

icall  =  0 

93 

if  (  typeC  i-t-1)  .ne.0  )  return 

99 

if  (  iscmd(symb( i+1) .command, ncmd) .gt .0  )  return 

100 

i  =  i  +  1 

101 

i n i t prm  =  symb ( i ) 

102 

lermin  =  3 

103 

return 

104 

c 

105 

c 

end 

command 

106 

c 

107 

180 

cont inue 

108 

lermin  =  I 

109 

relurn 

1 10 

c 

1 1 1 

c 

al  ib 

command 

1 12 

c 

113 

200 

coni inue 

114 

call  zmovechr  (  buf,  0,  libuf,  0,  nl ib 

1  15 

neul  ib  =  0 

116 

205 

coni inue 

1 17 

if  (  type ( i  +  1 ). gt .  1  )  go  lo  210 

1 18 

i  =  i  +  1 

119 

call  symove  (  libuf,  neul ib,  symbli' 

120 

if  (  type( i+1) .eq. 1  )  go  to  205 

121 

call  zmovechr  (  libuf,  neul  ib, 

122 

neul ib  =  neul ib  +  1 

123 

go  lo  205 

124 

210 

cont inue 

125 

call  zmovechr  (  libuf,  neul ib,  buf. 

126 

nl  ib  =  nl ib  +  neul ib 

127 

go  to  20 

128 

c 

129 

c 

savef 

command 

130 

c 

131 

220 

cont inue 

132 

if  (  type( i+1) .ne.0  )  go  to  225 

133 

i  =  i  +  1 

134 

call  savef  (  symbCi),  funct,  nl  ) 

135 

go  to  20 

136 

225 

cont inue 

137 

call  rnsgtor  (  "file  name  required". 

138 

go  to  20 

139 

c 

140 

c 

saves 

command 

141 

c 

142 

230 

coni inue 

143 

if  (  type( i+1 ) .ne . 0  )  go  lo  225 

144 

i  »  i  +  1 

145 

call  savef  (  symb(i),  prog,  npl  ) 

146 

go  to  20 

147 

c 

148 

c 

xeq  command 

149 

c 

150 

240 

coni inue 

151 

lermin  =  2 

152 

if  (  lype ( i+1 ) . gt . 0  )  return 

153 

if  (  iscmdtsymbt i+1) .command, ncmd) .gt . 

154 

i  =  i  +  1 

155 

ini tprm  =  symb ( i ) 

156 

termin  *  4 

157 

return 

158 

c 

159 

c 

dl  ib 

command 

160 

c 

161 

260 

coni inue 

162 

if  (  lype ( i+1 ). ne . 0  )  go  to  20 

1  ) 

nl  ib  ) 

) 


return 
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163 

164 

165 

166 
16? 
163 
169 
1?0 
171 
1?2 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 


i  =  i  +  1 

call  remove  (  symb(i).  li^uf.  nl  ib  ) 
go  to  260 
c 

c  kb  in  command 
c 

280  continue 

keep  =  1 
go  to  20 
c 

c  sub  command 
c 

300  continue 

ical 1  =  1 

if  (  type ( i+1 ) . ne . 0  )  go  to  305 

call  readfct  (  ious.  symb(i+l).  prog.  npl. 
if  (  ires. It. 0  )  go  to  305 
i  »  i  +  1 
go  to  20 

305  cont inue 

call  gelfct  (  prog,  npl  ) 
go  to  20 

c  symbol  with  more  than  8  characters 
c 

4B0  coni inue 

call  msqtor  C  "symbol  too  long".  15  1 
go  to  20 
c 

c  set  number  of  solutions 
c 

420  cont i nue 

n  =  symb(i) 
go  to  20 
c 

end 


ires 
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1  subroutine  shiftmsg  (  line.  1 en  ) 

2  implicit  integer(a-z) 

3  integer  line(l).  len 

4  c 

5  c  shiftmsg  shifts  off  the  first  symbol  from  line.  len  is  the 

6  c  number  of  characters  in  line. 

7  c 

8  parameter  (  blank=40b  ) 

9  c 

10  sbK  =  zscanlne  <  line.  0.  len.  blank  ) 

11  loc  =  zscanleq  (  line.  sbk.  len-sbk.  blank  ) 

12  nbk  =  zscanlne  (  line,  sbk+loc.  len-sbk-loc.  blank  ) 

13  len  =  len  -  sbk  -  loc  -  nbk 

14  call  zmovechr  C  line.  0.  line,  sbk+loc+nbk.  len  ) 

15  return 

16  c 

1?  end 
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subroutine  symove  (  destin.  loc*  source  ) 

implicit  inleger(a-z) 

integer  destin(l).  loc*  source 


1 

2 

3 


4 

c 

5 

c 

move  characters  from  source  to  destin(loc) 

unt 

6 

c 

detected.  loc  is  the  offset  in  characters 

of 

7 

c 

a 

parameter  (  blank=4Qb  ) 

9 

parameter  C  cpu=8  ) 

10 

c 

l  i 

n  =  zscanleq  (  source*  0*  cpw.  blank 

) 

12 

call  zmovechr  (  dest  in*  loc.  source. 

0* 

13 

loc  =  loc  +  n 

14 

return 

15 

c 

16 

end 

i  1  a  blank  is 
dest  in 


n  ) 
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subroutine  driver 


1 

2  c 


3 

c 

the  driver 

of  the  solver  generated  program 

4 

c 

5 

external 

f 

6 

complex 

f,  y (10241 

7 

integer 

1(30),  id C 10) -  inlen,  mapset,  plotmode. 

plotid#  of 

8 

integer 

grids,  dbwg.  first 

9 

complex 

c(30) .  p( 1024, 100) 

10 

real 

x( 1024) ,  u( 1024) 

11 

common 

/parmbk/  1.  c,  p,  y 

12 

common 

/plotbk/  id.  mapset.  plotmode,  plot  id. 

x#  u* 

13 

. 

xsizo,  ysize.  inlen.  grids,  nx 

:#  nu 

14 

common 

/miscbk/  n,  maxit.  knoun.  h.  epl.  ep2. 

of 

15 

namel ist 

/initbk/  n,  maxit,  known,  h.  epl.  ep2> 

:y 

% 

Q- 

O 

o 

16 

• 

mapset.  plotmode,  plotid.  x.  u.  xsize.  ysize 

1? 

. 

id.  grids,  dbug.  nx.  nu 

18 

parameter  (  npmax=1024.  nrmax=25  ) 

19 

real 

root re (npmax, nr max) .  root a i (npmax, nr max) 

20 

c 

2 , 

c 

def aul t  val ues 

22 

c 

23 

data 

maxit  /100/ 

24 

data 

knoun  /0/ 

25 

data 

h  /.5/ 

26 

data 

epl  /0 .  / 

27 

data 

ep2  /0 .  / 

28 

data 

of  /0/ 

29 

data 

mapset  /0/ 

30 

data 

plot mode  /0/ 

31 

data 

plotid  /0/ 

32 

data 

xsize  /6.0/ 

33 

data 

ysize  /4.5/ 

34 

data 

id  /'"solver"/ 

35 

data 

grids  /0/ 

36 

data 

dbug  /0/ 

37 

data 

nx  /npmax/ 

38 

data 

nu  /nrmax/ 

39 

c 

40 

integer 

line (10).  symb(l).  type(l) 

41 

integer 

narnef  (2) 

42 

data 

namef  /  5r+root.  4rrout  / 

43 

c 

44 

c 

initial izat ions 

45 

c 

46 

call  setfile  (  namef,  2  ) 

47 

call  dropfile  (  namef (1)  ) 

48 

cal  1  msg 

flag  (  msg,  2  ) 

49 

c 

50 

c 

read  input. 

i  f  any 

51 

c 

52 

i f  (  msg 

. eq . 0  )  go  to  15 

53 

cal  1 

msgfrr  (  line.  len.  80  ) 

54 

cal  1 

getnumb  (  symb.  type,  ns,  line.  len.  1  ) 

-50- 


55 

if  (  symb( 1) .ne. "tty"  >  go  to  20 

56 

call  mprompt  (  1  ) 

57 

read  (  59.  initbk  ) 

53 

call  mprompt  <  0,  0  ) 

59 

go  to  15 

60 

20 

cont inue 

61 

call  open  (  5.  symb(l).  0.  inlen  ) 

62 

if  (  inlen. It. .1  )  go  to  35 

63 

read  (  5.  initbK  ) 

64 

c 

65 

c 

create  and  write  to  output  file,  if  of=l  or  2 

66 

c 

67 

15 

cont inue 

68 

if  (  of.lt. 1  .or.  of.gt.2  1  go  to  35 

69 

name f (2)  =  400brt'namef (2)  +  140b 

70 

10 

cont inue 

71 

namef (2)  =  namef (2)  +  1 

72 

call  open  (  1.  namef (2).  0.  len  ) 

73 

cal  1  close  (  1  ) 

74 

if  (  len.gt.0  1  go  to  10 

75 

call  create  (  6,  namef (2).  2.  -1  ) 

76 

write  (  59.  103  )  namef (2) 

77 

write  (  6.  104  )  n.  maxit.  Known,  h.  epl.  ep2.  (y( 

70 

c 

79 

c 

checK  information 

80 

c 

31 

35 

con  t i nue 

82 

n  =  min0  (  n.  1024  ) 

83 

if  (  of.gt.2  .and.  pi otmode . eq .0  )  plotmode  =  1 

84 

iplot  =  pi otmode 

85 

if  (  iplot. gt. 3  .or.  mapset.gt.3  )  n  =  1 

86 

np  =  nx 

87 

nr  =  nw 

88 

if  (  np*nr . le. 25600  .and.  mapset.gt.3  )  go  to  40 

89 

np  =  npmax 

90 

nr  =  nr  max 

91 

c 

92 

c 

main  loop 

93 

c 

94 

40 

cont inue 

95 

call  wsetup 

96 

if  (  p 1 otmode . gt . 0  .and.  first. eq. 1  ) 

97 

.  call  rootplot  (  rootre.  rootai.  np.  nr  ) 

98 

if  (  1(301.11.0  )  go  to  70 

99 

f  irst  »  1 

100 

if  (  iplot. gt. 3  )  go  to  60 

101 

num  =  n 

102 

call  muller  (  Known,  num.  y.  maxit.  h»  epl.  ep2. 

103 

if  (  of.gt.2  )  go  to  40 

104 

if  (  of.eq.2  )  go  to  50 

105 

write  (  59.  101  )  (  y(i),  i  =  1.  num+Known  ) 

106 

if  (  kount . gt . max i t  )  write  (  59.  102  )  Kount 

107 

if  (  of.eq.0  )  go  to  40 

108 

50 

cont inue 

i  ■ 1 . n) 


kount  ) 


-51- 


103  write  (  6-101)  (y(i).  i=l.  num+knoun  ) 

110  if  (  Kount . gl . max i t  )  urile  C  6.  102)  Kounl 

111  go  to  40 
1  12  c 

113  c  evaluating  function  instead  of  finding  roots 

114  c 

115  60  continue 

1 16  ixind  =  ixind  +  i 

117  y(l)  =  f  C  x(ixind)  ) 

118  go  to  40 

1 19  c 

120  c  terminal ing 

121  c 

122  70  cont inue 

123  if  (  iplot.le.0  )  call  exit  (  dbug  ) 

124  plot mode  =  iplot 

125  call  rootplot  (  rootre.  rootai.  np.  nr  ) 

126  call  donepl 

127  call  exit  (  dbug  ) 

128  c 

129  101  format  (  "roots:''  ✓  <  *  (  ''.el9.ll."  .  ■.el9.ll.1’  )*  )  ) 

130  102  format  (  “does  not  converge  after  ”.  i5.  “  iterations"  ) 

131  103  format  <  a8.  "  is  the  output  file"  ) 

132  104  format  (  11  n  i4  /  “max it  »*.  i4  /  "Knoun  14/ 

133  .  “h  =".  f 15 . 7  /  “ep 1  =  ".  fl5.7  /  "ep2  fl5.7  ✓ 

134  .  “y  =".  2(2f 15.7.  ".  ")  /  (  7x.  2C2fl5.7.  ".•))) 

135  end 
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1  subroutine  curves  (  x,  y,  np,  mark,  nc  ) 

2  rool  x(  1) ,  y t 1024,  1 ) 

3  integer  np,  mark,  nc 

4  c 

5  c  draw  curves 

6  c 

?  do  10  i  =  1 ,  nc 

call  curve  (  x(l),  y(l,i),  np,  mark  ) 

9  10  cont inue 

10  return 

1 1  end 
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subroutine  destin  (  plotmode,  plot  id  ) 
integer  plotmode,  plot  id 


1 

2 


3 

c 

4 

c 

plot  mode  ■  1  ■=>  fr£0 

5 

c 

=  2  =>  pr tpl t 

6 

c 

=  3  =>  IK 

7 

c 

pi ot id  for  <  f r80> 

=  >  id 

8 

c 

for  <prtpl t> 

=  > • 1  ine 

width 

9 

c 

for  <tK> 

■*>  baud 

rate 

10 

c 

1 1 

mode  =  plolmode 

12 

i f  (  mode . 1 1  .  1  ) 

return 

13 

if  (  mode.gt.3  ) 

mode  = 

mod  ( 

14 

go  to  (  10,  20. 

30  ) »  mode 

15 

c 

16 

c 

f r80  file 

1? 

c 

18 

10  continue 

19 

call  Keep30  ( 

11  roots". 

3  ) 

20 

cal  1  f r80 id  ( 

plot  id. 

1,  1 

21 

call  pits 

22 

return 

23 

c 

24 

c 

pr  inter 

25 

c 

26 

20  con t inue 

27 

call  prtplt  ( 

plot  id  ) 

28 

return 

29 

c 

30 

c 

t  eK  tron ix 

31 

c 

32 

30  cont inue 

33 

cal  1  tK  (  0.  . 

l*pl ot  id 

.  59 

34 

return 

35 

c 

36 

end 

mode,  4  )  +  1 
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1 

2 

3 

4 

5  c 
S  c 
7  c 

a  c 

9  c 

10  c 

1 1  c 

12 

13 

14 

15  c 

16 
17 
13 

19 

20 
21 
22 

23 

24 

25 

26  c 
2  7 

28 

29 

30  c 

31 

32 

33 

34  c 

35 

36 
3? 

38  c 

39 

40 

41 

42  c 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54  c 

55 


subroutine  maptype  (  mapset.  npoint,  xmin.  xmax.  ymin,  ymax. 

.  xsize.  ysize  ) 

integer  mapset.  npoint 

real  xmin.  xmax,  ymin,  ymax,  xsize.  ysize 

mapsel  =0  =>  1  inear-1  inear 

=  1  =>  lineai — log 

=2  =>  log- linear 

=3  ->  log-log 

=4  =>  contour 

parameter  C  nusp  = 16384  ) 
common  uorK(nusp) 

data  scale  /"scale11/ 

map  =  mapset 

i f  (  map . 1 t . 0  )  map  =  0 
if  (  map.gt.4  )  map  =  0 
nxsize  =  (  xsize+1  )  /  2 
nysize  =  (  ysize+1  )  /  2 

xinc  =  stepsz  (  xmin,  xmax,  2.*nxsize,  xd  ) 

yinc  =  stepsz  (  ymin,  ymax,  2.>Knysize.  yd  ) 

xcyc  =  xsize  /  ncyc  (  xmin,  xmax,  xdec  ) 

ycyc  =  ysize  /  ncyc  (  ymin,  ymax,  ydec  ) 

go  to  (  10,  20,  30.  40,  50  ) .  map+1 

10  cont inue 

call  graf  (  xmin-xd.  scale,  xmax+xd.  ymin-yd.  scale,  ymox+yd  ) 
return 

20  cont inue 

call  ylog  (  xmin-xd,  xinc.  ydec.  ycyc  1 
return 

30  con  t  inue 

cell  xlog  (  xdoc,  xcyc,  ymin-yd.  yinc  1 
return 

40  cont inue 

call  loglog  (  xdec,  xcyc.  ydec.  ycyc  1 
return 

50  cont  inue 

xinc  =  xinc  -  2.  *  xd  /  maxB  (  2*nxsize.4  ) 

yinc  =  yinc  -  2.  *  xd  /  max0  (  2*nysize,4  ) 

call  graf  (  xmin,  xinc.  xmax.  ymin,  yinc.  ymax  ) 
call  bcomon  (  nusp  ) 

call  coni  in  (  0,  "solid",  "labels",  1.  10  ) 

call  coni  in  (  1.  "dash",  “nolabels".  1,  8  ) 

call  coni  in  (  2.  "chndot",  "nolabels".  1.  5  ) 

call  coni  in  (  3.  “dot",  "nolabels".  1,  3  ) 
call  conang  (60.  ) 
return 

end 
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1 

2 

3 

4  c 

5  c 

6  c 

7  c 

8  c 

9  c 

10  c 

1 1  c 

12  c 

13  c 

14  c 

15  c 
IS  c 

1 7  c 

18  c 

19  c 

20  c 

21  c 

22  c 

23  c 

24 

25 

26  c 

27  c 

28 

29 

30 

31 

32  c 

33 

34 

35  c 

36  c 

37 

38 

39 

40 

41 

42 

43 

44 

45 
4b 

47 

48 

49 

50 

51 

52  c 

53 

54 


subroutine  mullor  (Kn.  n.  rts.  maxit.  step.  epl.  ep2.  fn.  Kount) 
complex  r ts ( 1 ) 
complex  fn 

this  subroutine  uses  the  muller's  method  to  find  the  roots  of 

a  funct ion 

parameters: 

Kn  -  number  of  roots  previously  computed,  normally  zero 
n  -  input:  number  of  roots  desired 

output:  number  of  roots  found 

rts  -  an  array  of  initial  estimates  of  all  desired  roots 
set  to  zero  if  no  better  estimates  are  available 
max  it  -  maximum  number  of  iterations  per  root  allowed 
step  -  initial  distance  between  x’s 
epl  -  relative  error  tolerance  on  x(i) 
ep2  -  error  tolerance  on  f(x(i)) 

fn  -  fn(z)  is  a  external  complex  function  which,  for  given 

z.  returns  f (z) 

Kount  -  number  of  iterations  of  the  last  root.  if  Kount 
is  equal  to  maxit.  the  routine  will  stop  at  once 

complex  rt.  h.  delfpr.  frtdef.  lambda,  delf.  dfprlm.  num.  den 
complex  g.  sqr.  frt.  frtprv 

ini t ial  izal ion 

epsl  -  amaxl  (epl.  l.e-12) 
eps2  =  amaxl  (ep2.  l.e-20) 
ibeg  =  Kn  +  1 
iend  =  Kn  +  n 

do  100  i  3  ibeg.  iend 
Kount  =  0 

compute  first  three  estimates  for  root  as 

rts(i)+slep.  rts(i)-step.  and  rts(i) 

1  h  =  step 

rl  =rts(i)  +h 
index  =  1 
go  to  70 

10  del f  pr  =  f  r tdef 

rt  =  rts(  i)  -  h 
index  =  2 
go  to  70 

20  frtprv  =  frtdef 

delfpr  =  frtprv  -  delfpr 

r t  =  rts(  i ) 
index  =  3 
go  to  70 

30  index  =  4 

1 ambda  -  1 . 

compute  next  estimate  for  root 

40  delf  =  frtdef  -  frtprv 

dfprlm  =  delfpr  *  lambda 
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55 

56 

5 7 

58 

59 

60 
61 
62 

63 

64 

65 

66 
6? 
68 

69 

70 

71 

72 

73 

74 

75 
7  6 
77 
79 

79 

80 
81 
82 

83 

84 

85 

86 
8? 
88 

89 

90 

91 

92 

93 

94 

95 

96 
9? 


num  =  -frldef  *  (1.  +  lambda)  *  2. 
g  =  (1.  +  lambda  *  2.)  *  del f  -  lambda  *  dfprlm 
sqr  =  g  *  g  +  2.  *  num  *  lambda  *  (delf  -  dfprlm) 
sqr  =  csqrt (sqr) 
den  =  g  +  sqr 

if  (real (g)  *  real (sqr)  +  aimag(g)  *  aimag(sqr)  .It.  8.) 
1  den  =  g  -  sqr 

if  (cabs (den)  . 1 epl)  den  =  1. 

1  arnbda  =  num  /  den 
frtprv  =  frldef 
del f  pr  =  del f 
h  =  h  *  1 arnbda 
rl  =  rt  +  h 

if  (Kount  .gl.  maxil)  go  to  200 
c 

70  Kount  =  Kount  +  1 
frl  =  fn  (rt) 
frtdef  =  frt 

do  7 1  j  =  1 ,  i  - 1 
den  =  rt  -  rts(j) 

if  (cabs(den)  .It.  eps2)  go  to  79 
frtdef  =  frtdef  /  den 

7 1  cont inue 

75  go  to  (10,  20.  30.  83),  index 

79  r t s ( i )  =  rt  +  .001 
go  to  1 

c 

c  checK  for  convergence 

80  if  (cabs(h)  .It.  epsl  *  cabs(rt))  go  to  103 

if  (amax 1 (cabs ( f r t ) ,  cabs ( f rldef ) )  .11.  eps2)  go  to  100 
c 

c  checK  for  divergence 

if  (cabs ( f rtdef )  .It.  10.  *  cabs(frtprv) )  go  to  40 
h  =  .5  *  h 
lambda  =  .5  *  lambda 
rt  =  rt  -  h 
go  to  70 
100  rts(i)  =  rt 
return 

200  cont  .nue 
n  =  i 
return 

end 
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1  integer  function  ncyc  (  vmin.  vmax.  vdec  ) 

2  real  vmin.  vmax.  vdec 

3  c 

4  c  return  the  number  of  cycles  needed  for  log  scale 

5  c 

6  parameter  (  one =. 999999999999  ) 

7  parameter  (  maxcyc=10  ) 

8  parameter  (  base=l.&-20  ) 

9  c 

10  delta  *  0. 

11  if  (  'v'min.lt.  1.  )  delta  =  1. 

12  min  =  alogl0  (  amax 1 (vmin . base)  )  -  delta 

13  max  =  alog!0  (  amax 1 ( vmax. base)  )  +  one 

14  ncyc  =  max0  (  max-min.  1  ) 

15  vdec  =  10.  **  min 

16  if  (  vdec.gt .base*10.  .or.  ncyc . 1 t . maxcyc  )  return 

17  ncyc  =  rnaxcyc 

18  vdec  =  10.  **  (  max-maxcyc  ) 

19  return 

20  c 

2 1  end 
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1 

subroutine  rootplct  C  rootre 

.  rootai.  npneu.  nrneu  ) 

2 

real 

rootre (npneu. nrneu) , 

rootai (npneu. nrneu) 

3 

c 

4 

c 

plots  root  re  vs.  x  and/or  u  and 

rootai  vs.  x  and/or  u 

5 

c 

rootre  and 

rootai  are  the  real  and  imaginary  part  of  the 

6 

c 

roots  (or 

the  value  of  the  function,  if  so.  stored  in 

7 

c 

rootre ( 1 ) 

and  rootai (1)).  npneu 

is  the  maximum  number  of 

8 

c 

points  and 

nrneu  is  the  maximum 

number  of  roots. 

9 

c 

10 

integer 

1(30).  id (10).  inlen.  mapset.  plot  mode,  plot  id. 

1 1 

integer 

grids 

12 

real 

x ( 1024) .  u( 1024) 

13 

comp! ex 

y( 1024) 

14 

complex 

c (30) ,  p(1024, 100) 

15 

common 

/parmbK/  1.  c.  p. 

y 

16 

common 

/plotbK/  id.  mapset.  plotmode,  plot  id,  x.  u. 

17 

. 

xsize.  ysize.  inlen,  grids 

18 

common 

/miscbK/  n.  maxit. 

Knoun,  h,  ep 1 .  ep2.  of 

19 

c 

20 

real 

psize(2) 

21 

parameter  (  b  ig=  10  .  **  108,  smal 1 =- 13 . **100  ) 

22 

integer 

first,  npoint.  lino(ll) 

23 

data 

remin.  aimin.  xmin. 

umin  /4*big/ 

24 

data 

remax,  aimax,  xr.iax. 

umax  /4*small/ 

25 

data 

1  ine (11)  /lh$/ 

26 

data 

ps ize  /1 1 .  .8. 5/ 

27 

da  t  a 

nx  70/.  nu  /25630/ 

28 

c 

29 

integer 

labfre(2).  labfai(2).  labtre(2).  labtai(2) 

30 

integer 

labure(2).  labuai(2).  labx(2) 

31 

integer 

labpu(2).  labax(2) . 

1 abau(2> 

32 

data 

labfre  /“real.(  f ) " . 

II  u  , 

33 

data 

labfai  /"imcg(f)" . 

II  n  y 

34 

data 

labpu  /“paramote " , 

■r  u“/ 

35 

data 

labax  /  argument". 

“  X*/ 

36 

data 

labau  /"argument". 

"  u”/ 

37 

data 

1 abx  /"paramete". 

nr  x"/ 

38 

data 

labure  /“real(roo". 

•t)  */ 

39 

data 

1 abua i  /" imag(roo“ , 

*1)  */ 

40 

data 

labtre  /"  11  "/ 

41 

data 

labtai  /"  ",  ■  "/ 

42 

c 

43 

c 

initial izat ions 

44 

c 

45 

if  (  first. ns. 0  )  go  to  5 

46 

f irst  = 

1 

47 

nr  *  min0  <  n+Knoun*  nrneu  ) 

48 

xsz  =  xsize 

49 

ysz  =  ysize 

50 

c 

51 

c 

check  if  page  size  should  be  8.5 

"xl  1  ”  or  1 1  "x8.5" 

52 

c 

53 

if  (  xs 

z.eq.ysz  )  go  to  5 

54 

if  (  xs 

z.lt.ysz  .and.  ysz.le. 

psize(2)-1.5  )  go  to  5 
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55 

if 

(  xsz.ge.ysz  .and.  xsz . gt .ps ize (2) -  1 . 5  )  go  to  5 

56 

tenip  *  ps  ize(  1 ) 

57 

ps  ize( 1 )  3  ps  ize(2) 

58 

psize(2)  3  ternp 

59 

c 

60 

c 

store 

the  roots  of  the  function  value,  y 

61 

c 

62 

5 

cont inue 

63 

if 

(  1(30). It. 0  )  first  =  first  +  1 

64 

if 

(  firsl.gt.2  )  go  to  20 

65 

if 

(  mapset.gt.3  )  go  to  12 

66 

nx 

•  nx  +  1 

67 

if 

(  nx.gt.npneu  )  return 

68 

do 

10  i  =  1 .  nr 

69 

rootre(nx.i)  =  real  (  y(i)  ) 

70 

rootai(nx.i)  3  aimag  (  y(i)  ) 

71 

if  (  rernin . gt . rootre (nx.  i )  )  remin  3  rootre(nx.i) 

72 

if  (  remax. 1 t ,rootre(nx. i)  )  remax  3  rootre(nx.i) 

73 

if  (  aimin.gt .rootai (nx. i)  )  aimin  3  rootai(nx.i) 

74 

if  (  a imax . 1 l . roota i (nx. i )  )  aimax  3  roolai(nx.i) 

75 

10 

coni inue 

76 

return 

77 

c 

78 

12 

cont  inue 

79 

if  (  nu.  1 1  .nrneu  )  go  to  14 

80 

nx  3  nx  +  1 

81 

nu  3  0 

82 

14 

cont inue 

83 

nu  =  nu  +  1 

84 

rootre(nx.nu)  3  real  (  y(l)  ) 

85 

rootai (nx.nu)  3  aimag  (  y(l)  ) 

86 

if  (  remin. gt .rootre(nx.nu)  )  remin  3  roplre(nx.nu) 

87 

if  (  rernax.  1 1  .rootre(nx.nu)  )  remax  3  rootre  (nx.  nu) 

88 

if  (  aimin.gt .rootai (nx.nu)  )  aimin  3  rootai (nx.nu) 

89 

if  (  aimax. 1 t .rootai (nx.nu)  )  aimax  3  rootai (nx.nu) 

90 

return 

91 

c 

92 

c 

f  ind 

the  minimum  and  maximum  values  of  x  and  u 

93 

c 

94 

20 

cont inue 

95 

npoint  =  min0  (  nx.  npneu  ) 

96 

do  30  j  •  1.  r point 

97 

if  (  xmin.gt.x(j)  )  xmin  3  x(j) 

98 

if  (  xmax.lt.x(j)  )  xmax  3  x( j) 

99 

30 

cont inue 

100 

if  (  mapset.lt.4  )  go  to  40 

101 

do  35  j  3  1 .  nu 

102 

if  (  umin.gt.u(j)  )  umin  3  u(j) 

103 

if  (  umax.ll.u(j)  )  umax  3  u(j) 

104 

35 

cont inue 

105 

if  (  remax-remin. 1 t . 1 .e-25  )  none  3  1 

106 

if  (  a imax-a imin . 1 t . 1 . e-25  )  noai  3  1 

107 

remin  3  umin 

103 

aimin  3  umin 
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remax  =  umax 
a i max  *>  umax 


figure  out  the  title  and  name  for  the  axes 
40  cont inue 

if  (  plotmode. 1 t .4  .and.  mapset.lt. 4  ) 
if  <  pi o t mode  1 1 t . 4  )  go  to  50 


1  17 

cal  1 

move 

( 

labx.  labax.  2  ) 

1  18 

cal  1 

move 

c 

1 abure. 

labfre.  2  ) 

1 19 

cal  1 

move 

c 

labuai. 

1 abf a i .  2  ) 

120 

if  ( 

mapse  t 

1  t  .4  ) 

go  to  60 

121 

cal  1 

move 

( 

1 abure. 

1 abau.  2  ) 

122 

cal  1 

move 

t 

1 abua i . 

labau.  2  ) 

123 

cal  1 

move 

c 

1 ab  tre. 

labfre.  2  ) 

124 

cal  1 

move 

c 

labtai. 

1 abf  a i .  2  ) 

125 

go  to  60 

126 

50 

cont inue 

127 

cal  1 

move 

c 

1 abtre. 

1 abure.  2  ) 

128 

cal  1 

move 

< 

labtai. 

1 abua i .  2  ) 

129 

cal  1 

move 

( 

1 abure. 

labpu.  2  ) 

130 

cal  1 

move 

c 

1 abua i . 

labpu.  2  ) 

go  to  60 


getting  ready  to  plot 

60  cont inue 

call  setbox  (  id.  80  ) 
call  gfsize  (  3.  2000000b  ) 
call  destin  (  plot  mode.,  plot  id  ) 
cal  1  nobrdr 


plot  the  input  file,  if  there  is  one 

if  (  inlen.le.0  )  go  to  90 

call  page  (  psize(l).  psize(2)  ) 
call  tablet  (  "Itset".  “long"  ) 
rewind  5 
70  cont inue 

read(5.1010)  (  line(i).  i = 1 , 10  ) 

if  (  ieof(5).ne.0  )  go  to  80 
call  1 1 1  ine  (  line  ) 
go  to  70 
80  cont inue 

call  endtab  (  0  ) 
cal  1  dendpl  (  0  ) 

98  cont inue 

xsz  *  aminl  (  xsz.  psize(l)-1.5  ) 
ysr  «  aminl  (  ysz.  psize(2)-1.5  ) 
call  page  (  psize(l).  psize(2)  ) 

plot  the  real (root)  vs.  x 

xr.or  «  (  ps  ize(  1 )  *. 5-xsz  )*.75 
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163 

if  C  xsor.lt. .75  )  xsor  =  (  psize(l)- 

xsz  )*.5 

1 6*4 

ysor  ■  (  psire(2)*.5-ysz  )*.5 

165 

if  (  ysor.lt. .5  )  ysor  ■  (  ps ize (2) -ysz  )*.5 

166 

call  physor  (  xsor.  ysor  ) 

167 

c 

168 

coll  title  (  labtre,  16.  labx.  16.  labuire.  16.  xsz.  ysz  ) 

169 

call  maptype  (  mapset.  npoint.  xmin.  xmax.  remin.  remax. 

170 

. 

ksz.  ysz  ) 

171 

call  d frame 

172 

if  (  grids. ne.0  )  call  grid  (  -grids. 

-grids  ) 

173 

if  (  mapset. It. 4  )  go  to  100 

174 

if  (  nore.eq.l  )  go  to  120 

175 

cal  1  conmaK  (  rootre.  npneu.  nrneu. 

"scale"  ) 

176 

call  contur  (  4.  “labels",  "drau"  ) 

177 

go  to  110 

173 

100 

cont  inue 

179 

call  curves  (  x.  rootre.  npoint.  -1 

.  nr  ) 

180 

1 10 

cont inue 

181 

call  endgr  (  0  ) 

182 

c 

183 

c 

plot  imag i nary (root)  vs.  x 

184 

c 

185 

120 

cont  inue 

186 

xsor  =  psize(l)*.5  +  (  psizef  l)*.5-xsz 

)*.5 

187 

if  (  xsor . 1 l . . 5+ps ize ( 1 ) *. 5  )  xsor  = 

(  psize(l)-xsz  )*.5 

188 

ysor  =  ps ize (2) *. 5  +  (  psize(2)*.5-ysz 

)*.5 

189 

if  (ysor . 1 t . .5+psize(2)*.5  )  ysor  »  ( 

psize(2)-ysz  )*.5 

190 

call  physor  (  xsor.  ysor  ) 

191 

if  (  xsor . 1 t . ps ize( 1 ) *. 5  .and.  ysor.lt 

.psize(2)*.5  ) 

192 

call  dendpl  (  0  ) 

193 

c 

194 

call  title  (  lablai.  16,  labx.  16,  labuai.  16.  xsz,  ysz  1 

195 

call  maptype  (  mapset.  npoint.  xmin.  xmax.  aimin.  aimax. 

196 

xsz.  ysz  ) 

197 

call  dframe 

198 

if  (  grids. ne.0  )  call  grid  (  -grids. 

-grids  ) 

199 

if  (  mapset. It. 4  1  go  to  130 

200 

if  (  noai.eq.l  )  go  to  140 

201 

call  conmaK  (  rootai.  npneu,  nrneu. 

"scale"  ) 

202 

call  contur  (  4.  “labels",  “drau"  ) 

203 

go  to  140 

204 

130 

cont inue 

205 

call  curves  (  x.  rootai,  npoint.  -1 

,  nr  ) 

206 

140 

cont inue 

207 

call  endgr  (  0  ) 

208 

call  dendpl  (  0  ) 

209 

c 

210 

return 

211 

c 

212 

1010 

format  (  10a8  ) 

213 

c 

214 

end 
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1  subroutine  set  file  (names,  n) 

2  c 

3  c  subroutine  setfilo  appends  the  current  suffix  (channel)  to 

4  c  all  the  given  files. 

5  c 

6  c  parameters: 

7  c  names  -  an 

8  c  in 

9c  n  -  an 

10  c 

11  c  1  ibrary  used : 

12  c  f  or  1 1  ib 

13  c 

14  c 

15  integer  names(n).  n 

16  integer  userno.  account,  df.  suffix,  i 

1?  c 

18  c  set  unit  59  to  be  the  controller  or  terminal 

19  c 

20  call  msgl inK  (59.  1) 

21  c 

22  c  gel  the  suf  f  ix 

23  c 

24  call  userinfo  (userno.  account,  df.  suffix) 

25  c 

26  c  change  the  1 ef l- just i f ied  suffix  to  r igh t - just i f ied 

27  c 

28  suffix  =  suffix  /  4000000000000000000b 

29  c 

30  c  append  the  suffix  to  al 1  file  names 

31  c 

32  do  10  i  =  1.  n 

33  names(i)  =  400b  *  names(i)  +  suffix 

34  10  continue 

35  c 

36  return 

3?  end 


one-dimensional  integer  array  with  file  names 
it.  all’ names  are  limited  to  7  characters, 
integer  that  indicates  the  number  of  files. 
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uDOD^icr>aix^wro 


1  reol  function  stepsz  (  vmin.  umax,  size,  delta  ) 

real  vmin.  vmax.  size,  delta 

c 

c  find  step  size  for  a  graph 
c 

delta  =  amaxl  (  0.001.  (vmax-vmin)*.05  ) 
stepsz  =  (  vmax-vmin+2 . *del ta  )  /  max0  (  int(size).  4 
return 
end 
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fNiro^Tin^N-aocn 


subroutine  wsetup 

complex  c (30) ,  p( 1024, 108),  y(1024) 

integer  1 (30) 

common  /parrnbK/  1.  c,  p.  y 
data  1 (30)  /{/ 

1 (30)  =  1 t30)  -  1 

return 

end 
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